1 /* *****************************************************************
2 MESQUITE -- The Mesh Quality Improvement Toolkit
3
4 Copyright 2004 Sandia Corporation and Argonne National
5 Laboratory. Under the terms of Contract DE-AC04-94AL85000
6 with Sandia Corporation, the U.S. Government retains certain
7 rights in this software.
8
9 This library is free software; you can redistribute it and/or
10 modify it under the terms of the GNU Lesser General Public
11 License as published by the Free Software Foundation; either
12 version 2.1 of the License, or (at your option) any later version.
13
14 This library is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 Lesser General Public License for more details.
18
19 You should have received a copy of the GNU Lesser General Public License
20 (lgpl.txt) along with this library; if not, write to the Free Software
21 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22
23 diachin2@llnl.gov, djmelan@sandia.gov, mbrewer@sandia.gov,
24 pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov
25
26 ***************************************************************** */
27 // -*- Mode : c++; tab-width: 3; c-tab-always-indent: t; indent-tabs-mode: nil; c-basic-offset: 3 -*-
28 //
29 // SUMMARY:
30 // USAGE:
31 //
32 // AUTHOR: Thomas Leurent <tleurent@mcs.anl.gov>
33 // ORG: Argonne National Laboratory
34 // E-MAIL: tleurent@mcs.anl.gov
35 //
36 // ORIG-DATE: 12-Nov-02 at 18:05:56
37 // LAST-MOD: 7-May-03 at 13:28:47 by Thomas Leurent
38 //
39 // DESCRIPTION:
40 // ============
41 /*! \file MsqMeshEntityTest.cpp
42
43 Unit testing of various functions in the MsqMeshEntity class.
44
45 \author Michael Brewer
46 \author Thomas Leurent
47
48 */
49 // DESCRIP-END.
50 //
51
52
53 #include "MsqMeshEntity.hpp"
54 #include "Vector3D.hpp"
55 #include "PatchData.hpp"
56 #include "PatchDataInstances.hpp"
57 #include <math.h>
58 #include <iostream>
59 #include <sstream>
60 #include "UnitUtil.hpp"
61
62 using namespace MBMesquite;
63 using std::cout;
64 using std::endl;
65
66 class MsqMeshEntityTest : public CppUnit::TestFixture
67 {
68
69 private:
70 CPPUNIT_TEST_SUITE(MsqMeshEntityTest);
71 CPPUNIT_TEST (test_hex_vertices);
72 CPPUNIT_TEST (test_centroid_tri);
73 CPPUNIT_TEST (test_centroid_quad);
74 CPPUNIT_TEST (test_centroid_hex);
75 CPPUNIT_TEST (test_unsigned_area);
76 CPPUNIT_TEST (test_unsigned_area_poly);
77 CPPUNIT_TEST (test_unsigned_area_tet);
78 CPPUNIT_TEST (test_unsigned_area_pyr);
79 CPPUNIT_TEST (test_unsigned_area_pri);
80 CPPUNIT_TEST (test_unsigned_area_hex);
81 CPPUNIT_TEST (test_all_nodes);
82 CPPUNIT_TEST (test_check_element_orientation_linear);
83 CPPUNIT_TEST (test_check_element_orientation_quadratic);
84 CPPUNIT_TEST_SUITE_END();
85
86 void test_all_nodes( EntityTopology type, unsigned num_nodes );
87
88 private:
89 PatchData oneHexPatch;
90 PatchData oneTetPatch;
91 PatchData oneQuadPatch;
92 PatchData oneTriPatch;
93 Vector3D e1, e2, e3;
94 double tolEps;
95 public:
setUp()96 void setUp()
97 {
98 tolEps=1.e-12;
99
100 // sets up the unit vectors
101 e1.set(1,0,0);
102 e2.set(0,1,0);
103 e3.set(0,0,1);
104
105 MsqPrintError err(cout);
106
107 // creates empty Patch
108 create_one_hex_patch(oneHexPatch, err); CPPUNIT_ASSERT(!err);
109 create_one_tet_patch(oneTetPatch, err); CPPUNIT_ASSERT(!err);
110 create_one_tri_patch(oneTriPatch, err); CPPUNIT_ASSERT(!err);
111 create_one_quad_patch(oneQuadPatch, err); CPPUNIT_ASSERT(!err);
112 }
113
tearDown()114 void tearDown()
115 {
116 destroy_patch_with_domain( oneTriPatch );
117 destroy_patch_with_domain( oneQuadPatch );
118 }
119
120 public:
MsqMeshEntityTest()121 MsqMeshEntityTest()
122 { }
123
test_hex_vertices()124 void test_hex_vertices()
125 {
126 MsqPrintError err(cout);
127 // prints out the vertices.
128 const MsqVertex* ideal_vertices = oneHexPatch.get_vertex_array(err); CPPUNIT_ASSERT(!err);
129 size_t num_vtx = oneHexPatch.num_nodes();
130 CPPUNIT_ASSERT_EQUAL(size_t(8), num_vtx);
131
132 MsqVertex vtx;
133
134 vtx.set(1,1,1);
135 CPPUNIT_ASSERT_EQUAL(vtx, ideal_vertices[0]);
136
137 vtx.set(2,2,2);
138 CPPUNIT_ASSERT_EQUAL(vtx, ideal_vertices[6]);
139
140 vtx.set(1,2,2);
141 CPPUNIT_ASSERT_EQUAL(vtx, ideal_vertices[7]);
142 }
143
144 //! test the centroid of the first element in the Patch
test_centroid(PatchData & pd,Vector3D & correct)145 void test_centroid(PatchData &pd, Vector3D &correct)
146 {
147 MsqPrintError err(cout);
148 double eps = 1e-6;
149 Vector3D centroid;
150
151 MsqMeshEntity* elem = pd.get_element_array(err); CPPUNIT_ASSERT(!err);
152 elem->get_centroid(centroid, pd, err); CPPUNIT_ASSERT(!err);
153
154 // cout << "centroid: "<< centroid <<endl;
155 // cout << "correct: "<< correct <<endl;
156
157 for (int i=0; i<3; ++i)
158 CPPUNIT_ASSERT_DOUBLES_EQUAL( centroid[i], correct[i], eps);
159 }
160
test_centroid_tri()161 void test_centroid_tri()
162 {
163 Vector3D correct(1.5, 1+1/(2.0*sqrt(3.0)), 1.0);
164 test_centroid(oneTriPatch, correct);
165 }
166
test_centroid_quad()167 void test_centroid_quad()
168 {
169 Vector3D correct(1.5, 1.5, 1.0);
170 test_centroid(oneQuadPatch, correct);
171 }
172
test_centroid_hex()173 void test_centroid_hex()
174 {
175 Vector3D correct(1.5, 1.5, 1.5);
176 test_centroid(oneHexPatch, correct);
177 }
178
179
test_unsigned_area()180 void test_unsigned_area()
181 {
182 MsqPrintError err(cout);
183 MsqMeshEntity* tri = oneTriPatch.get_element_array(err);
184 CPPUNIT_ASSERT(!err);
185 CPPUNIT_ASSERT(fabs(tri->compute_unsigned_area(oneTriPatch,err)
186 -(sqrt(3.0)/4.0)) < tolEps);
187 MsqMeshEntity* quad = oneQuadPatch.get_element_array(err);
188 CPPUNIT_ASSERT(!err);
189 CPPUNIT_ASSERT(fabs(quad->compute_unsigned_area(oneQuadPatch,err)
190 -1.0) < tolEps);
191 }
192
193 void test_unsigned_area_poly();
194 void test_unsigned_area_tet();
195 void test_unsigned_area_pyr();
196 void test_unsigned_area_pri();
197 void test_unsigned_area_hex();
198 void test_all_nodes();
199
200 void test_unsigned_area_common( EntityTopology type,
201 const double* coords,
202 double expected );
203
204 void test_check_element_orientation_linear();
205 void test_check_element_orientation_quadratic();
206 void test_check_element_orientation( EntityTopology type, int nodes );
207 };
208
209
210 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION(MsqMeshEntityTest, "MsqMeshEntityTest");
211 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION(MsqMeshEntityTest, "Unit");
212
213 const size_t conn[] = { 0, 1, 2, 3, 4, 5, 6, 7 };
214 const bool fixed[] = { false, false, false, false, false, false, false, false };
215
test_unsigned_area_poly()216 void MsqMeshEntityTest::test_unsigned_area_poly()
217 {
218 const double coords[] = { 0, 0, 0,
219 1, 0, 0,
220 1, 1, 0,
221 0.5, 1.5, 0,
222 0, 1, 0 };
223 size_t n_vtx = 5;
224 EntityTopology type = POLYGON;
225 MsqError err;
226
227 PatchData pd;
228 pd.fill( n_vtx, coords, 1, &type, &n_vtx, conn, fixed, err );
229 ASSERT_NO_ERROR(err);
230
231 double a = pd.element_by_index(0).compute_unsigned_area( pd, err );
232 ASSERT_NO_ERROR(err);
233 CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.25, a, 1e-8 );
234 }
235
test_unsigned_area_common(EntityTopology type,const double * coords,double expected)236 void MsqMeshEntityTest::test_unsigned_area_common( EntityTopology type,
237 const double* coords,
238 double expected )
239 {
240 MsqError err;
241 PatchData pd;
242
243 pd.fill( TopologyInfo::corners( type ), coords, 1, type, conn, fixed, err );
244 ASSERT_NO_ERROR(err);
245
246 double a = pd.element_by_index(0).compute_unsigned_area( pd, err );
247 ASSERT_NO_ERROR(err);
248
249 CPPUNIT_ASSERT_DOUBLES_EQUAL( expected, a, 1e-8 );
250 }
251
252
test_unsigned_area_tet()253 void MsqMeshEntityTest::test_unsigned_area_tet()
254 {
255 const double coords[] = { 0, 0, 0,
256 1, 0, 0,
257 0, 1, 0,
258 0, 0, 1 };
259 test_unsigned_area_common( TETRAHEDRON, coords, 1.0/6.0 );
260 }
261
test_unsigned_area_pyr()262 void MsqMeshEntityTest::test_unsigned_area_pyr()
263 {
264 const double coords[] = { 0, 0, 0,
265 1, 0, 0,
266 1, 1, 0,
267 0, 1, 0,
268 0, 0, 1 };
269 test_unsigned_area_common( PYRAMID, coords, 1.0/3.0 );
270
271 const double pyr_coords[] = {-1, -1, -1,
272 1, -1, -1,
273 1, 1, -1,
274 -1, 1, -1,
275 0, 0, 0 };
276 test_unsigned_area_common( PYRAMID, pyr_coords, 4.0/3.0 );
277 }
278
test_unsigned_area_pri()279 void MsqMeshEntityTest::test_unsigned_area_pri()
280 {
281 const double coords[] = { 0, 0, 0,
282 1, 0, 0,
283 0, 1, 0,
284 0, 0, 1,
285 1, 0, 1,
286 0, 1, 1 };
287 test_unsigned_area_common( PRISM, coords, 0.5 );
288
289 const double tet_coords[] = { 0, 0, 0,
290 1, 0, 0,
291 0, 1, 0,
292 0, 0, 1,
293 0, 0, 1,
294 0, 0, 1 };
295 test_unsigned_area_common( PRISM, tet_coords, 1.0/6.0 );
296 }
297
test_unsigned_area_hex()298 void MsqMeshEntityTest::test_unsigned_area_hex()
299 {
300 const double coords[] = { 0, 0, 0,
301 1, 0, 0,
302 1, 1, 0,
303 0, 1, 0,
304 0, 0, 1,
305 1, 0, 1,
306 1, 1, 1,
307 0, 1, 1 };
308 test_unsigned_area_common( HEXAHEDRON, coords, 1.0 );
309
310 const double coords2[] = { 0, 0, 0,
311 2, 0, 0,
312 2, 2, 0,
313 0, 2, 0,
314 0, 0, 2,
315 2, 0, 2,
316 2, 2, 2,
317 0, 2, 2 };
318 test_unsigned_area_common( HEXAHEDRON, coords2, 8.0 );
319
320 const double pyr_coords[] = {-1,-1, 0,
321 1,-1, 0,
322 1, 1, 0,
323 -1, 1, 0,
324 0, 0, 1,
325 0, 0, 1,
326 0, 0, 1,
327 0, 0, 1 };
328 test_unsigned_area_common( HEXAHEDRON, pyr_coords, 4.0/3.0 );
329 }
330
test_all_nodes(EntityTopology type,unsigned num_nodes)331 void MsqMeshEntityTest::test_all_nodes( EntityTopology type, unsigned num_nodes )
332 {
333 const unsigned num_vtx = 27;
334 double coords[3*num_vtx] = {0.0};
335 size_t conn[num_vtx];
336 for (size_t i = 0; i < num_vtx; ++i)
337 conn[i] = i;
338 bool fixed[num_vtx] = {false};
339 CPPUNIT_ASSERT(num_nodes <= num_vtx);
340
341 MsqError err;
342 PatchData pd;
343 size_t n = num_nodes;
344 pd.fill( num_nodes, coords, 1, &type, &n, conn, fixed, err );
345 ASSERT_NO_ERROR(err);
346
347 MsqMeshEntity& elem = pd.element_by_index(0);
348 NodeSet all = elem.all_nodes(err);
349 ASSERT_NO_ERROR(err);
350 CPPUNIT_ASSERT_EQUAL( num_nodes, all.num_nodes() );
351 CPPUNIT_ASSERT( all.have_any_corner_node() );
352 bool mid_edge, mid_face, mid_reg;
353 TopologyInfo::higher_order( type, num_nodes, mid_edge, mid_face, mid_reg, err );
354 ASSERT_NO_ERROR(err);
355 CPPUNIT_ASSERT_EQUAL( mid_edge, !!all.have_any_mid_edge_node() );
356 CPPUNIT_ASSERT_EQUAL( mid_face, !!all.have_any_mid_face_node() );
357 CPPUNIT_ASSERT_EQUAL( mid_reg, !!all.have_any_mid_region_node() );
358
359 }
360
test_all_nodes()361 void MsqMeshEntityTest::test_all_nodes()
362 {
363 test_all_nodes( TRIANGLE, 3 );
364 test_all_nodes( TRIANGLE, 4 );
365 test_all_nodes( TRIANGLE, 6 );
366 test_all_nodes( TRIANGLE, 7 );
367
368 test_all_nodes( QUADRILATERAL, 4 );
369 test_all_nodes( QUADRILATERAL, 5 );
370 test_all_nodes( QUADRILATERAL, 8 );
371 test_all_nodes( QUADRILATERAL, 9 );
372
373 test_all_nodes( TETRAHEDRON, 4 );
374 test_all_nodes( TETRAHEDRON, 5 );
375 test_all_nodes( TETRAHEDRON, 10 );
376 test_all_nodes( TETRAHEDRON, 11 );
377 test_all_nodes( TETRAHEDRON, 8 );
378 test_all_nodes( TETRAHEDRON, 9 );
379 test_all_nodes( TETRAHEDRON, 14 );
380 test_all_nodes( TETRAHEDRON, 15 );
381
382 test_all_nodes( HEXAHEDRON, 8 );
383 test_all_nodes( HEXAHEDRON, 9 );
384 test_all_nodes( HEXAHEDRON, 20 );
385 test_all_nodes( HEXAHEDRON, 21 );
386 test_all_nodes( HEXAHEDRON, 14 );
387 test_all_nodes( HEXAHEDRON, 15 );
388 test_all_nodes( HEXAHEDRON, 26 );
389 test_all_nodes( HEXAHEDRON, 27 );
390 }
391
test_check_element_orientation_linear()392 void MsqMeshEntityTest::test_check_element_orientation_linear()
393 {
394 const EntityTopology types[] = { TRIANGLE,
395 QUADRILATERAL,
396 TETRAHEDRON,
397 PYRAMID,
398 PRISM,
399 HEXAHEDRON };
400 const int num_types = sizeof(types)/sizeof(types[0]);
401
402 for (int i = 0; i < num_types; ++i) {
403 test_check_element_orientation( types[i], TopologyInfo::corners(types[i]) );
404 }
405 }
406
test_check_element_orientation_quadratic()407 void MsqMeshEntityTest::test_check_element_orientation_quadratic()
408 {
409 struct ElemType {
410 EntityTopology topo;
411 unsigned nodes;
412 };
413
414 const ElemType types[] = {
415 { TRIANGLE, 6 },
416 { QUADRILATERAL, 8 },
417 { QUADRILATERAL, 9 },
418 { TETRAHEDRON, 10 } };
419 const int num_types = sizeof(types)/sizeof(types[0]);
420
421 for (int i = 0; i < num_types; ++i) {
422 test_check_element_orientation( types[i].topo, types[i].nodes );
423 }
424 }
425
test_check_element_orientation(EntityTopology type,int nodes)426 void MsqMeshEntityTest::test_check_element_orientation( EntityTopology type,
427 int nodes )
428 {
429 // get an ideal element
430 MsqError err;
431 PatchData pd;
432 create_ideal_element_patch( pd, type, nodes, err );
433 ASSERT_NO_ERROR(err);
434 CPPUNIT_ASSERT_EQUAL( (size_t)1, pd.num_elements() );
435 CPPUNIT_ASSERT_EQUAL( (size_t)nodes, pd.num_nodes() );
436 MsqMeshEntity& elem = pd.element_by_index(0);
437 CPPUNIT_ASSERT_EQUAL( (size_t)nodes, elem.node_count() );
438 CPPUNIT_ASSERT_EQUAL( type, elem.get_element_type() );
439 const size_t* conn = elem.get_vertex_index_array();
440
441 // test that ideal element is not reported as inverted
442 int inverted, tested;
443 elem.check_element_orientation( pd, inverted, tested, err );
444 ASSERT_NO_ERROR(err);
445 CPPUNIT_ASSERT_EQUAL( 0, inverted );
446 CPPUNIT_ASSERT( tested > 0 );
447
448 bool mids[4] = {false};
449 TopologyInfo::higher_order( type, nodes, mids[1], mids[2], mids[3], err );
450 MSQ_ERRRTN(err);
451
452 // invert element at each vertex and test
453 Vector3D centroid;
454 elem.get_centroid( centroid, pd, err );
455 ASSERT_NO_ERROR(err);
456 for (int i = 0; i < nodes; ++i) {
457 unsigned dim, num;
458 TopologyInfo::side_from_higher_order( type, nodes, i, dim, num, err );
459 ASSERT_NO_ERROR(err);
460 const Vector3D old_pos = pd.vertex_by_index( conn[i] );
461 Vector3D new_pos = old_pos;
462 if (dim == TopologyInfo::dimension(type)) {
463 // move mid-element node 3/4 of the way to corner 0
464 new_pos += 3*pd.vertex_by_index( conn[0] );
465 new_pos *= 0.25;
466 }
467 else if (dim == 0) { // if a corner vertex
468 if (type == TRIANGLE || type == TETRAHEDRON) {
469 // move tri/tet vertex past opposite side of element
470 new_pos += 2*(centroid - old_pos);
471 }
472 else if (mids[1]) {
473 // if have mid-edge nodes move 3/4 of the way to center vertex
474 new_pos += 3*centroid;
475 new_pos *= 0.25;
476 }
477 else {
478 // move vertex past centroid
479 new_pos += 1.5*(centroid - old_pos);
480 }
481 }
482 else {
483 // otherwise move vertex past centroid
484 new_pos += 2.5*(centroid - old_pos);
485 }
486
487 pd.set_vertex_coordinates( new_pos, conn[i], err );
488 ASSERT_NO_ERROR(err);
489
490 // test that element is inverted
491 inverted = tested = 0;
492 elem.check_element_orientation( pd, inverted, tested, err );
493 ASSERT_NO_ERROR(err);
494 std::ostringstream str;
495 str << TopologyInfo::short_name(type) << nodes
496 << " Vertex " << i
497 << " (Dimension " << dim
498 << " Index " << num << ")";
499 CppUnit::Message m( "MsqMeshEntity failed to detect inverted element" );
500 m.addDetail( str.str() );
501 ASSERT_MESSAGE( m, inverted > 0 );
502
503 // move vertex back to ideal position
504 pd.set_vertex_coordinates( old_pos, conn[i], err );
505 ASSERT_NO_ERROR(err);
506 }
507 }
508