00001 /************************************************************************************* 00002 * MechSys - A C++ library to simulate (Continuum) Mechanical Systems * 00003 * Copyright (C) 2005 Dorival de Moraes Pedroso <dorival.pedroso at gmail.com> * 00004 * Copyright (C) 2005 Raul Dario Durand Farfan <raul.durand at gmail.com> * 00005 * * 00006 * This file is part of MechSys. * 00007 * * 00008 * MechSys is free software; you can redistribute it and/or modify it under the * 00009 * terms of the GNU General Public License as published by the Free Software * 00010 * Foundation; either version 2 of the License, or (at your option) any later * 00011 * version. * 00012 * * 00013 * MechSys is distributed in the hope that it will be useful, but WITHOUT ANY * 00014 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A * 00015 * PARTICULAR PURPOSE. See the GNU General Public License for more details. * 00016 * * 00017 * You should have received a copy of the GNU General Public License along with * 00018 * MechSys; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, * 00019 * Fifth Floor, Boston, MA 02110-1301, USA * 00020 *************************************************************************************/ 00021 00022 #ifndef MECHSYS_CONREC_H 00023 #define MECHSYS_CONREC_H 00024 00025 #include <stdio.h> 00026 #include <math.h> 00027 #include <vector> 00028 00029 int conrec(double **d, 00030 int ilb, 00031 int iub, 00032 int jlb, 00033 int jub, 00034 double *x, 00035 double *y, 00036 int nc, 00037 double const *z, 00038 std::vector<double> & X, 00039 std::vector<double> & Y); 00040 00041 #define xsect(p1,p2) (h[p2]*xh[p1]-h[p1]*xh[p2])/(h[p2]-h[p1]) 00042 #define ysect(p1,p2) (h[p2]*yh[p1]-h[p1]*yh[p2])/(h[p2]-h[p1]) 00043 #define min(x,y) (x<y?x:y) 00044 #define max(x,y) (x>y?x:y) 00045 00046 /* 00047 Copyright (c) 1996-1997 Nicholas Yue 00048 00049 This software is copyrighted by Nicholas Yue. This code is base on the work of 00050 Paul D. Bourke CONREC.F routine 00051 00052 The authors hereby grant permission to use, copy, and distribute this 00053 software and its documentation for any purpose, provided that existing 00054 copyright notices are retained in all copies and that this notice is included 00055 verbatim in any distributions. Additionally, the authors grant permission to 00056 modify this software and its documentation for any purpose, provided that 00057 such modifications are not distributed without the explicit consent of the 00058 authors and that existing copyright notices are retained in all copies. Some 00059 of the algorithms implemented by this software are patented, observe all 00060 applicable patent law. 00061 00062 IN NO EVENT SHALL THE AUTHORS OR DISTRIBUTORS BE LIABLE TO ANY PARTY FOR 00063 DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT 00064 OF THE USE OF THIS SOFTWARE, ITS DOCUMENTATION, OR ANY DERIVATIVES THEREOF, 00065 EVEN IF THE AUTHORS HAVE BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00066 00067 THE AUTHORS AND DISTRIBUTORS SPECIFICALLY DISCLAIM ANY WARRANTIES, INCLUDING, 00068 BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A 00069 PARTICULAR PURPOSE, AND NON-INFRINGEMENT. THIS SOFTWARE IS PROVIDED ON AN 00070 "AS IS" BASIS, AND THE AUTHORS AND DISTRIBUTORS HAVE NO OBLIGATION TO PROVIDE 00071 MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. 00072 */ 00073 00074 //============================================================================= 00075 // 00076 // CONREC is a contouring subroutine for rectangularily spaced data. 00077 // 00078 // It emits calls to a line drawing subroutine supplied by the user 00079 // which draws a contour map corresponding to real*4data on a randomly 00080 // spaced rectangular grid. The coordinates emitted are in the same 00081 // units given in the x() and y() arrays. 00082 // 00083 // Any number of contour levels may be specified but they must be 00084 // in order of increasing value. 00085 // 00086 // As this code is ported from FORTRAN-77, please be very careful of the 00087 // various indices like ilb,iub,jlb and jub, remeber that C/C++ indices 00088 // starts from zero (0) 00089 // 00090 //============================================================================= 00091 int conrec(double **d, 00092 int ilb, 00093 int iub, 00094 int jlb, 00095 int jub, 00096 double *x, 00097 double *y, 00098 int nc, 00099 double const *z, 00100 std::vector<double> & X, 00101 std::vector<double> & Y) 00102 // d ! matrix of data to contour 00103 // ilb,iub,jlb,jub ! index bounds of data matrix 00104 // x ! data matrix column coordinates 00105 // y ! data matrix row coordinates 00106 // nc ! number of contour levels 00107 // z ! contour levels in increasing order 00108 { 00109 int m1,m2,m3,case_value; 00110 double dmin,dmax; 00111 double x1=0,x2=0,y1=0,y2=0; 00112 register int i,j,k,m; 00113 double h[5]; 00114 int sh[5]; 00115 double xh[5],yh[5]; 00116 //=========================================================================== 00117 // The indexing of im and jm should be noted as it has to start from zero 00118 // unlike the fortran counter part 00119 //=========================================================================== 00120 int im[4] = {0,1,1,0},jm[4]={0,0,1,1}; 00121 //=========================================================================== 00122 // Note that castab is arranged differently from the FORTRAN code because 00123 // Fortran and C/C++ arrays are transposed of each other, in this case 00124 // it is more tricky as castab is in 3 dimension 00125 //=========================================================================== 00126 int castab[3][3][3] = 00127 { 00128 { 00129 {0,0,8},{0,2,5},{7,6,9} 00130 }, 00131 { 00132 {0,3,4},{1,3,1},{4,3,0} 00133 }, 00134 { 00135 {9,6,7},{5,2,0},{8,0,0} 00136 } 00137 }; 00138 for (j=(jub-1);j>=jlb;j--) { 00139 for (i=ilb;i<=iub-1;i++) { 00140 double temp1,temp2; 00141 temp1 = min(d[i][j],d[i][j+1]); 00142 temp2 = min(d[i+1][j],d[i+1][j+1]); 00143 dmin = min(temp1,temp2); 00144 temp1 = max(d[i][j],d[i][j+1]); 00145 temp2 = max(d[i+1][j],d[i+1][j+1]); 00146 dmax = max(temp1,temp2); 00147 if (dmax>=z[0]&&dmin<=z[nc-1]) { 00148 for (k=0;k<nc;k++) { 00149 if (z[k]>=dmin&&z[k]<=dmax) { 00150 for (m=4;m>=0;m--) { 00151 if (m>0) { 00152 //============================================================= 00153 // The indexing of im and jm should be noted as it has to 00154 // start from zero 00155 //============================================================= 00156 h[m] = d[i+im[m-1]][j+jm[m-1]]-z[k]; 00157 xh[m] = x[i+im[m-1]]; 00158 yh[m] = y[j+jm[m-1]]; 00159 } else { 00160 h[0] = 0.25*(h[1]+h[2]+h[3]+h[4]); 00161 xh[0]=0.5*(x[i]+x[i+1]); 00162 yh[0]=0.5*(y[j]+y[j+1]); 00163 } 00164 if (h[m]>0.0) { 00165 sh[m] = 1; 00166 } else if (h[m]<0.0) { 00167 sh[m] = -1; 00168 } else 00169 sh[m] = 0; 00170 } 00171 //================================================================= 00172 // 00173 // Note: at this stage the relative heights of the corners and the 00174 // centre are in the h array, and the corresponding coordinates are 00175 // in the xh and yh arrays. The centre of the box is indexed by 0 00176 // and the 4 corners by 1 to 4 as shown below. 00177 // Each triangle is then indexed by the parameter m, and the 3 00178 // vertices of each triangle are indexed by parameters m1,m2,and 00179 // m3. 00180 // It is assumed that the centre of the box is always vertex 2 00181 // though this isimportant only when all 3 vertices lie exactly on 00182 // the same contour level, in which case only the side of the box 00183 // is drawn. 00184 // 00185 // 00186 // vertex 4 +-------------------+ vertex 3 00187 // | \ / | 00188 // | \ m-3 / | 00189 // | \ / | 00190 // | \ / | 00191 // | m=2 X m=2 | the centre is vertex 0 00192 // | / \ | 00193 // | / \ | 00194 // | / m=1 \ | 00195 // | / \ | 00196 // vertex 1 +-------------------+ vertex 2 00197 // 00198 // 00199 // 00200 // Scan each triangle in the box 00201 // 00202 //================================================================= 00203 for (m=1;m<=4;m++) { 00204 m1 = m; 00205 m2 = 0; 00206 if (m!=4) 00207 m3 = m+1; 00208 else 00209 m3 = 1; 00210 case_value = castab[sh[m1]+1][sh[m2]+1][sh[m3]+1]; 00211 if (case_value!=0) { 00212 switch (case_value) { 00213 //=========================================================== 00214 // Case 1 - Line between vertices 1 and 2 00215 //=========================================================== 00216 case 1: 00217 x1=xh[m1]; 00218 y1=yh[m1]; 00219 x2=xh[m2]; 00220 y2=yh[m2]; 00221 break; 00222 //=========================================================== 00223 // Case 2 - Line between vertices 2 and 3 00224 //=========================================================== 00225 case 2: 00226 x1=xh[m2]; 00227 y1=yh[m2]; 00228 x2=xh[m3]; 00229 y2=yh[m3]; 00230 break; 00231 //=========================================================== 00232 // Case 3 - Line between vertices 3 and 1 00233 //=========================================================== 00234 case 3: 00235 x1=xh[m3]; 00236 y1=yh[m3]; 00237 x2=xh[m1]; 00238 y2=yh[m1]; 00239 break; 00240 //=========================================================== 00241 // Case 4 - Line between vertex 1 and side 2-3 00242 //=========================================================== 00243 case 4: 00244 x1=xh[m1]; 00245 y1=yh[m1]; 00246 x2=xsect(m2,m3); 00247 y2=ysect(m2,m3); 00248 break; 00249 //=========================================================== 00250 // Case 5 - Line between vertex 2 and side 3-1 00251 //=========================================================== 00252 case 5: 00253 x1=xh[m2]; 00254 y1=yh[m2]; 00255 x2=xsect(m3,m1); 00256 y2=ysect(m3,m1); 00257 break; 00258 //=========================================================== 00259 // Case 6 - Line between vertex 3 and side 1-2 00260 //=========================================================== 00261 case 6: 00262 x1=xh[m3]; 00263 y1=yh[m3]; 00264 x2=xsect(m1,m2); 00265 y2=ysect(m1,m2); 00266 break; 00267 //=========================================================== 00268 // Case 7 - Line between sides 1-2 and 2-3 00269 //=========================================================== 00270 case 7: 00271 x1=xsect(m1,m2); 00272 y1=ysect(m1,m2); 00273 x2=xsect(m2,m3); 00274 y2=ysect(m2,m3); 00275 break; 00276 //=========================================================== 00277 // Case 8 - Line between sides 2-3 and 3-1 00278 //=========================================================== 00279 case 8: 00280 x1=xsect(m2,m3); 00281 y1=ysect(m2,m3); 00282 x2=xsect(m3,m1); 00283 y2=ysect(m3,m1); 00284 break; 00285 //=========================================================== 00286 // Case 9 - Line between sides 3-1 and 1-2 00287 //=========================================================== 00288 case 9: 00289 x1=xsect(m3,m1); 00290 y1=ysect(m3,m1); 00291 x2=xsect(m1,m2); 00292 y2=ysect(m1,m2); 00293 break; 00294 default: 00295 break; 00296 } 00297 //============================================================= 00298 // Put your processing code here and comment out the printf 00299 //============================================================= 00300 //printf("%f %f %f %f %f\n",x1,y1,x2,y2,z[k]); 00301 X.push_back(x1); X.push_back(x2); 00302 Y.push_back(y1); Y.push_back(y2); 00303 } 00304 } 00305 } 00306 } 00307 } 00308 } 00309 } 00310 return 0; 00311 } 00312 00313 #endif // MECHSYS_CONREC_H 00314 00315 // vim:fdm=marker
1.4.7