|
MeshKit
1.0
|
|
MeshKit
1.0
|
00001 #include "meshkit/Mesh.hpp" 00002 00003 using namespace Jaal; 00004 00005 void quality_values( const vector<double> &quality) 00006 { 00007 cout << " MinValue : " << *min_element( quality.begin(), quality.end() ) << endl; 00008 cout << " MaxValue : " << *max_element( quality.begin(), quality.end() ) << endl; 00009 cout << " Mean Value : " << quality[quality.size()/2] << endl; 00010 size_t nSize = quality.size(); 00011 double sum = 0.0; 00012 for( unsigned int i = 0; i < nSize; i++) 00013 sum += quality[i]; 00014 cout << " Avg Value : " << sum/(double)quality.size() << endl; 00015 } 00016 00017 void print_histogram( const vector<double> &quality, const string &header, int n ) 00018 { 00019 string filename; 00020 switch( n ) { 00021 case 3: 00022 filename ="./QData/T" + header + ".gnu"; 00023 break; 00024 case 4: 00025 filename = "./QData/Q" + header + ".gnu"; 00026 break; 00027 } 00028 00029 ofstream ofile( filename.c_str(), ios::out); 00030 if( ofile.fail() ) { 00031 cout << "Warning: Can not open file " << filename << endl; 00032 return; 00033 } 00034 00035 ofile << "set title \" QuadMesh Quality: " << header << "\"" << " font \"Helvetica,20\"" << endl; 00036 ofile << "set xlabel \"Quad Elements\" font \"Helvetica,15 \"" << endl; 00037 ofile << "set ylabel \"" << header << "\" font \"Helvetica,15\"" << endl; 00038 ofile << "set format y \"%11.3E\"" << endl; 00039 ofile << "set grid" << endl; 00040 ofile << "set terminal png" << endl; 00041 00042 switch( n ) { 00043 case 3: 00044 filename = "./QData/T" + header + ".data"; 00045 ofile << "set output \"" << "T" << header << ".png\"" << endl; 00046 ofile << "plot \"" << "T" << header << ".data\"" << " notitle with impulses" << endl; 00047 break; 00048 case 4: 00049 filename = "./QData/Q" + header + ".data"; 00050 ofile << "set output \"" << "Q" << header << ".png\"" << endl; 00051 ofile << "plot \"" << "Q" << header << ".data\"" << " notitle with impulses" << endl; 00052 break; 00053 } 00054 ofile.close(); 00055 00056 ofile.open( filename.c_str(), ios::out); 00057 if( ofile.fail() ) return; 00058 00059 for( unsigned int i = 0; i < quality.size(); i++) 00060 ofile << i << " " << quality[i] << endl; 00061 00062 } 00064 00065 void plot_all_quad_quality_measures( Jaal::Mesh *mesh ) 00066 { 00067 #ifdef HAVE_VERDICT 00068 VERDICT_REAL coords[4][3]; 00069 Point3D xyz; 00070 int num_nodes = 4; 00071 00072 size_t numfaces = mesh->getSize(2); 00073 00074 vector<double> quality( numfaces ); 00075 00077 00078 for( size_t i = 0; i < numfaces; i++) { 00079 Face *face = mesh->getFaceAt(i); 00080 for( int j = 0; j < 4; j++) { 00081 Vertex *vertex = face->getNodeAt(j); 00082 xyz = vertex->getXYZCoords(); 00083 coords[j][0] = xyz[0]; 00084 coords[j][1] = xyz[1]; 00085 coords[j][2] = xyz[2]; 00086 } 00087 quality[i] = v_quad_aspect( num_nodes, coords); 00088 } 00089 sort( quality.begin(), quality.end() ); 00090 print_histogram( quality, "AspectRatio", 4 ); 00091 cout << "Quality: AspectRatio " << endl; 00092 quality_values( quality ); 00093 00095 00096 for( size_t i = 0; i < numfaces; i++) { 00097 Face *face = mesh->getFaceAt(i); 00098 for( int j = 0; j < 4; j++) { 00099 Vertex *vertex = face->getNodeAt(j); 00100 xyz = vertex->getXYZCoords(); 00101 coords[j][0] = xyz[0]; 00102 coords[j][1] = xyz[1]; 00103 coords[j][2] = xyz[2]; 00104 } 00105 quality[i] = v_quad_skew( num_nodes, coords ); 00106 } 00107 sort( quality.begin(), quality.end() ); 00108 print_histogram( quality, "Skewness", 4 ); 00109 00111 00112 for( size_t i = 0; i < numfaces; i++) { 00113 Face *face = mesh->getFaceAt(i); 00114 for( int j = 0; j < 4; j++) { 00115 Vertex *vertex = face->getNodeAt(j); 00116 xyz = vertex->getXYZCoords(); 00117 coords[j][0] = xyz[0]; 00118 coords[j][1] = xyz[1]; 00119 coords[j][2] = xyz[2]; 00120 } 00121 quality[i] = v_quad_taper( num_nodes, coords ); 00122 } 00123 sort( quality.begin(), quality.end() ); 00124 print_histogram( quality, "Taper", 4 ); 00125 00127 00128 for( size_t i = 0; i < numfaces; i++) { 00129 Face *face = mesh->getFaceAt(i); 00130 for( int j = 0; j < 4; j++) { 00131 Vertex *vertex = face->getNodeAt(j); 00132 xyz = vertex->getXYZCoords(); 00133 coords[j][0] = xyz[0]; 00134 coords[j][1] = xyz[1]; 00135 coords[j][2] = xyz[2]; 00136 } 00137 quality[i] = v_quad_warpage( num_nodes, coords ); 00138 } 00139 sort( quality.begin(), quality.end() ); 00140 print_histogram( quality, "Warpage", 4 ); 00141 00142 00144 00145 for( size_t i = 0; i < numfaces; i++) { 00146 Face *face = mesh->getFaceAt(i); 00147 for( int j = 0; j < 4; j++) { 00148 Vertex *vertex = face->getNodeAt(j); 00149 xyz = vertex->getXYZCoords(); 00150 coords[j][0] = xyz[0]; 00151 coords[j][1] = xyz[1]; 00152 coords[j][2] = xyz[2]; 00153 } 00154 quality[i] = v_quad_area( num_nodes, coords ); 00155 } 00156 sort( quality.begin(), quality.end() ); 00157 print_histogram( quality, "Area", 4 ); 00158 cout << "Quality: Area " << endl; 00159 quality_values( quality ); 00160 00161 00163 00164 for( size_t i = 0; i < numfaces; i++) { 00165 Face *face = mesh->getFaceAt(i); 00166 for( int j = 0; j < 4; j++) { 00167 Vertex *vertex = face->getNodeAt(j); 00168 xyz = vertex->getXYZCoords(); 00169 coords[j][0] = xyz[0]; 00170 coords[j][1] = xyz[1]; 00171 coords[j][2] = xyz[2]; 00172 } 00173 quality[i] = v_quad_stretch( num_nodes, coords ); 00174 } 00175 sort( quality.begin(), quality.end() ); 00176 print_histogram( quality, "Stretch", 4 ); 00177 00179 00180 for( size_t i = 0; i < numfaces; i++) { 00181 Face *face = mesh->getFaceAt(i); 00182 for( int j = 0; j < 4; j++) { 00183 Vertex *vertex = face->getNodeAt(j); 00184 xyz = vertex->getXYZCoords(); 00185 coords[j][0] = xyz[0]; 00186 coords[j][1] = xyz[1]; 00187 coords[j][2] = xyz[2]; 00188 } 00189 quality[i] = v_quad_minimum_angle( num_nodes, coords ); 00190 } 00191 sort( quality.begin(), quality.end() ); 00192 print_histogram( quality, "MinAngle", 4 ); 00193 cout << "Quality: MinAngle (in Degrees) " << endl; 00194 quality_values( quality ); 00195 00197 00198 for( size_t i = 0; i < numfaces; i++) { 00199 Face *face = mesh->getFaceAt(i); 00200 for( int j = 0; j < 4; j++) { 00201 Vertex *vertex = face->getNodeAt(j); 00202 xyz = vertex->getXYZCoords(); 00203 coords[j][0] = xyz[0]; 00204 coords[j][1] = xyz[1]; 00205 coords[j][2] = xyz[2]; 00206 } 00207 quality[i] = v_quad_maximum_angle( num_nodes, coords ); 00208 } 00209 sort( quality.begin(), quality.end() ); 00210 print_histogram( quality, "MaxAngle", 4 ); 00211 cout << "Quality: MaxAngle (in Degrees) " << endl; 00212 quality_values( quality ); 00213 00215 00216 for( size_t i = 0; i < numfaces; i++) { 00217 Face *face = mesh->getFaceAt(i); 00218 for( int j = 0; j < 4; j++) { 00219 Vertex *vertex = face->getNodeAt(j); 00220 xyz = vertex->getXYZCoords(); 00221 coords[j][0] = xyz[0]; 00222 coords[j][1] = xyz[1]; 00223 coords[j][2] = xyz[2]; 00224 } 00225 quality[i] = v_quad_oddy( num_nodes, coords ); 00226 } 00227 sort( quality.begin(), quality.end() ); 00228 print_histogram( quality, "Oddy", 4 ); 00229 00231 00232 for( size_t i = 0; i < numfaces; i++) { 00233 Face *face = mesh->getFaceAt(i); 00234 for( int j = 0; j < 4; j++) { 00235 Vertex *vertex = face->getNodeAt(j); 00236 xyz = vertex->getXYZCoords(); 00237 coords[j][0] = xyz[0]; 00238 coords[j][1] = xyz[1]; 00239 coords[j][2] = xyz[2]; 00240 } 00241 quality[i] = v_quad_condition( num_nodes, coords ); 00242 } 00243 sort( quality.begin(), quality.end() ); 00244 print_histogram( quality, "ConditionNumber", 4 ); 00245 00247 00248 for( size_t i = 0; i < numfaces; i++) { 00249 Face *face = mesh->getFaceAt(i); 00250 for( int j = 0; j < 4; j++) { 00251 Vertex *vertex = face->getNodeAt(j); 00252 xyz = vertex->getXYZCoords(); 00253 coords[j][0] = xyz[0]; 00254 coords[j][1] = xyz[1]; 00255 coords[j][2] = xyz[2]; 00256 } 00257 quality[i] = v_quad_jacobian( num_nodes, coords ); 00258 } 00259 sort( quality.begin(), quality.end() ); 00260 print_histogram( quality, "Jacobian", 4 ); 00261 00263 00264 for( size_t i = 0; i < numfaces; i++) { 00265 Face *face = mesh->getFaceAt(i); 00266 for( int j = 0; j < 4; j++) { 00267 Vertex *vertex = face->getNodeAt(j); 00268 xyz = vertex->getXYZCoords(); 00269 coords[j][0] = xyz[0]; 00270 coords[j][1] = xyz[1]; 00271 coords[j][2] = xyz[2]; 00272 } 00273 quality[i] = v_quad_scaled_jacobian( num_nodes, coords ); 00274 } 00275 sort( quality.begin(), quality.end() ); 00276 print_histogram( quality, "ScaledJacobian", 4 ); 00277 00279 00280 for( size_t i = 0; i < numfaces; i++) { 00281 Face *face = mesh->getFaceAt(i); 00282 for( int j = 0; j < 4; j++) { 00283 Vertex *vertex = face->getNodeAt(j); 00284 xyz = vertex->getXYZCoords(); 00285 coords[j][0] = xyz[0]; 00286 coords[j][1] = xyz[1]; 00287 coords[j][2] = xyz[2]; 00288 } 00289 quality[i] = v_quad_shear( num_nodes, coords ); 00290 } 00291 sort( quality.begin(), quality.end() ); 00292 print_histogram( quality, "Shear", 4 ); 00293 00295 00296 for( size_t i = 0; i < numfaces; i++) { 00297 Face *face = mesh->getFaceAt(i); 00298 for( int j = 0; j < 4; j++) { 00299 Vertex *vertex = face->getNodeAt(j); 00300 xyz = vertex->getXYZCoords(); 00301 coords[j][0] = xyz[0]; 00302 coords[j][1] = xyz[1]; 00303 coords[j][2] = xyz[2]; 00304 } 00305 quality[i] = v_quad_shape( num_nodes, coords ); 00306 } 00307 sort( quality.begin(), quality.end() ); 00308 print_histogram( quality, "Shape", 4 ); 00309 00311 00312 for( size_t i = 0; i < numfaces; i++) { 00313 Face *face = mesh->getFaceAt(i); 00314 for( int j = 0; j < 4; j++) { 00315 Vertex *vertex = face->getNodeAt(j); 00316 xyz = vertex->getXYZCoords(); 00317 coords[j][0] = xyz[0]; 00318 coords[j][1] = xyz[1]; 00319 coords[j][2] = xyz[2]; 00320 } 00321 quality[i] = v_quad_relative_size_squared( num_nodes, coords ); 00322 } 00323 sort( quality.begin(), quality.end() ); 00324 print_histogram( quality, "RelativeSizeSquared", 4 ); 00325 00327 00328 for( size_t i = 0; i < numfaces; i++) { 00329 Face *face = mesh->getFaceAt(i); 00330 for( int j = 0; j < 4; j++) { 00331 Vertex *vertex = face->getNodeAt(j); 00332 xyz = vertex->getXYZCoords(); 00333 coords[j][0] = xyz[0]; 00334 coords[j][1] = xyz[1]; 00335 coords[j][2] = xyz[2]; 00336 } 00337 quality[i] = v_quad_shape_and_size( num_nodes, coords ); 00338 } 00339 sort( quality.begin(), quality.end() ); 00340 print_histogram( quality, "ShapeSize", 4 ); 00341 00343 00344 for( size_t i = 0; i < numfaces; i++) { 00345 Face *face = mesh->getFaceAt(i); 00346 for( int j = 0; j < 4; j++) { 00347 Vertex *vertex = face->getNodeAt(j); 00348 xyz = vertex->getXYZCoords(); 00349 coords[j][0] = xyz[0]; 00350 coords[j][1] = xyz[1]; 00351 coords[j][2] = xyz[2]; 00352 } 00353 quality[i] = v_quad_shear_and_size( num_nodes, coords ); 00354 } 00355 sort( quality.begin(), quality.end() ); 00356 print_histogram( quality, "ShearSize", 4 ); 00357 00359 00360 for( size_t i = 0; i < numfaces; i++) { 00361 Face *face = mesh->getFaceAt(i); 00362 for( int j = 0; j < 4; j++) { 00363 Vertex *vertex = face->getNodeAt(j); 00364 xyz = vertex->getXYZCoords(); 00365 coords[j][0] = xyz[0]; 00366 coords[j][1] = xyz[1]; 00367 coords[j][2] = xyz[2]; 00368 } 00369 quality[i] = v_quad_distortion( num_nodes, coords ); 00370 } 00371 sort( quality.begin(), quality.end() ); 00372 print_histogram( quality, "Distortion", 4 ); 00373 cout << "Quality: Distortion " << endl; 00374 quality_values( quality ); 00375 00376 00378 #endif 00379 } 00380 00381 void plot_all_tri_quality_measures( Jaal::Mesh *mesh ) 00382 { 00383 #ifdef HAVE_VERDICT 00384 VERDICT_REAL coords[3][3]; 00385 Point3D xyz; 00386 int num_nodes = 3; 00387 00388 size_t numfaces = mesh->getSize(2); 00389 00390 vector<double> quality( numfaces ); 00391 00393 00394 for( size_t i = 0; i < numfaces; i++) { 00395 Face *face = mesh->getFaceAt(i); 00396 for( int j = 0; j < num_nodes; j++) { 00397 Vertex *vertex = face->getNodeAt(j); 00398 xyz = vertex->getXYZCoords(); 00399 coords[j][0] = xyz[0]; 00400 coords[j][1] = xyz[1]; 00401 coords[j][2] = xyz[2]; 00402 } 00403 quality[i] = v_tri_aspect( num_nodes, coords); 00404 } 00405 sort( quality.begin(), quality.end() ); 00406 print_histogram( quality, "AspectRatio",3 ); 00407 00409 00410 for( size_t i = 0; i < numfaces; i++) { 00411 Face *face = mesh->getFaceAt(i); 00412 for( int j = 0; j < 3; j++) { 00413 Vertex *vertex = face->getNodeAt(j); 00414 xyz = vertex->getXYZCoords(); 00415 coords[j][0] = xyz[0]; 00416 coords[j][1] = xyz[1]; 00417 coords[j][2] = xyz[2]; 00418 } 00419 quality[i] = v_tri_area( num_nodes, coords ); 00420 } 00421 sort( quality.begin(), quality.end() ); 00422 print_histogram( quality, "Area" ,3); 00423 00425 00426 for( size_t i = 0; i < numfaces; i++) { 00427 Face *face = mesh->getFaceAt(i); 00428 for( int j = 0; j < 3; j++) { 00429 Vertex *vertex = face->getNodeAt(j); 00430 xyz = vertex->getXYZCoords(); 00431 coords[j][0] = xyz[0]; 00432 coords[j][1] = xyz[1]; 00433 coords[j][2] = xyz[2]; 00434 } 00435 quality[i] = v_tri_minimum_angle( num_nodes, coords ); 00436 } 00437 sort( quality.begin(), quality.end() ); 00438 print_histogram( quality, "MinAngle" ,3); 00439 00440 00442 00443 for( size_t i = 0; i < numfaces; i++) { 00444 Face *face = mesh->getFaceAt(i); 00445 for( int j = 0; j < 3; j++) { 00446 Vertex *vertex = face->getNodeAt(j); 00447 xyz = vertex->getXYZCoords(); 00448 coords[j][0] = xyz[0]; 00449 coords[j][1] = xyz[1]; 00450 coords[j][2] = xyz[2]; 00451 } 00452 quality[i] = v_tri_maximum_angle( num_nodes, coords ); 00453 } 00454 sort( quality.begin(), quality.end() ); 00455 print_histogram( quality, "MaxAngle" ,3); 00456 00458 00459 for( size_t i = 0; i < numfaces; i++) { 00460 Face *face = mesh->getFaceAt(i); 00461 for( int j = 0; j < 3; j++) { 00462 Vertex *vertex = face->getNodeAt(j); 00463 xyz = vertex->getXYZCoords(); 00464 coords[j][0] = xyz[0]; 00465 coords[j][1] = xyz[1]; 00466 coords[j][2] = xyz[2]; 00467 } 00468 quality[i] = v_tri_condition( num_nodes, coords ); 00469 } 00470 sort( quality.begin(), quality.end() ); 00471 print_histogram( quality, "ConditionNumber" ,3); 00472 00474 00475 for( size_t i = 0; i < numfaces; i++) { 00476 Face *face = mesh->getFaceAt(i); 00477 for( int j = 0; j < 3; j++) { 00478 Vertex *vertex = face->getNodeAt(j); 00479 xyz = vertex->getXYZCoords(); 00480 coords[j][0] = xyz[0]; 00481 coords[j][1] = xyz[1]; 00482 coords[j][2] = xyz[2]; 00483 } 00484 quality[i] = v_tri_scaled_jacobian( num_nodes, coords ); 00485 } 00486 sort( quality.begin(), quality.end() ); 00487 print_histogram( quality, "ScaledJacobian" ,3); 00488 00490 00491 for( size_t i = 0; i < numfaces; i++) { 00492 Face *face = mesh->getFaceAt(i); 00493 for( int j = 0; j < 3; j++) { 00494 Vertex *vertex = face->getNodeAt(j); 00495 xyz = vertex->getXYZCoords(); 00496 coords[j][0] = xyz[0]; 00497 coords[j][1] = xyz[1]; 00498 coords[j][2] = xyz[2]; 00499 } 00500 quality[i] = v_tri_shape( num_nodes, coords ); 00501 } 00502 sort( quality.begin(), quality.end() ); 00503 print_histogram( quality, "Shape" ,3); 00504 00506 00507 for( size_t i = 0; i < numfaces; i++) { 00508 Face *face = mesh->getFaceAt(i); 00509 for( int j = 0; j < 3; j++) { 00510 Vertex *vertex = face->getNodeAt(j); 00511 xyz = vertex->getXYZCoords(); 00512 coords[j][0] = xyz[0]; 00513 coords[j][1] = xyz[1]; 00514 coords[j][2] = xyz[2]; 00515 } 00516 quality[i] = v_tri_shape_and_size( num_nodes, coords ); 00517 } 00518 sort( quality.begin(), quality.end() ); 00519 print_histogram( quality, "ShapeSize" ,3); 00520 00522 00523 for( size_t i = 0; i < numfaces; i++) { 00524 Face *face = mesh->getFaceAt(i); 00525 for( int j = 0; j < 3; j++) { 00526 Vertex *vertex = face->getNodeAt(j); 00527 xyz = vertex->getXYZCoords(); 00528 coords[j][0] = xyz[0]; 00529 coords[j][1] = xyz[1]; 00530 coords[j][2] = xyz[2]; 00531 } 00532 quality[i] = v_tri_distortion( num_nodes, coords ); 00533 } 00534 sort( quality.begin(), quality.end() ); 00535 print_histogram( quality, "Distortion" ,3); 00536 00538 #endif 00539 }
1.7.6.1 
