Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
inks/MsqMeshEntity.cpp
Go to the documentation of this file.
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 //
28 // ORIG-DATE: 16-May-02 at 10:26:21
29 // LAST-MOD: 18-Jun-04 at 11:36:07 by Thomas Leurent
30 //
42 #include "Mesquite.hpp"
43 #include "MsqMeshEntity.hpp"
44 #include "MsqVertex.hpp"
45 #include "PatchData.hpp"
46 #include "MeshSet.hpp"
47 
48 #ifdef MSQ_USE_OLD_STD_HEADERS
49 # include <vector.h>
50 #else
51 # include <vector>
52  using std::vector;
53 #endif
54 
55 #ifdef MSQ_USE_OLD_IO_HEADERS
56 # include <ostream.h>
57 #else
58 # include <ostream>
59  using std::ostream;
60  using std::endl;
61 #endif
62 
63 namespace Mesquite {
64 
70 void Mesquite::MsqMeshEntity::get_vertex_indices(vector<size_t> &vertices) const
71 {
72  vertices.resize( vertex_count() );
73  memcpy( &vertices[0], vertexIndices, vertex_count() * sizeof(size_t) );
74 }
75 
83 void Mesquite::MsqMeshEntity::append_vertex_indices(vector<size_t> &vertex_list) const
84 {
85  vertex_list.insert(vertex_list.end(),
86  vertexIndices,
87  vertexIndices + vertex_count());
88 }
89 
90 void MsqMeshEntity::get_node_indices( vector<size_t>& indices ) const
91 {
92  indices.resize( node_count() );
93  memcpy( &indices[0], vertexIndices, node_count() * sizeof( size_t ) );
94 }
95 
96 void MsqMeshEntity::append_node_indices( vector<size_t>& indices ) const
97 {
98  indices.insert( indices.end(), vertexIndices, vertexIndices + node_count() );
99 }
100 
101 
106 void MsqMeshEntity::get_centroid(Vector3D &centroid, const PatchData &pd, MsqError &err) const
107 {
108  const MsqVertex* vtces = pd.get_vertex_array(err); MSQ_ERRRTN(err);
109  size_t nve = vertex_count();
110  for (size_t i=0; i<nve; ++i)
111  centroid += vtces[vertexIndices[i]];
112  centroid /= nve;
113 }
114 
115 
118 void Mesquite::MsqMeshEntity::compute_weighted_jacobian(PatchData &pd,
119  Vector3D &sample_point,
120  Vector3D jacobian_vectors[],
121  short &num_jacobian_vectors,
122  MsqError &err )
123 {
124  // v_v is just an alias for vertexIndices
125  const size_t* v_v = &vertexIndices[0];
126 
127  // vector<size_t> v_v;
128  // get_vertex_indices(v_v);
129  MsqVertex *vertices=pd.get_vertex_array(err); MSQ_ERRRTN(err);
130  switch (mType)
131  {
132  //Note:: For the linear tri case we do not use sample pt.
133  case TRIANGLE:
134  jacobian_vectors[0].set(vertices[v_v[1]]-vertices[v_v[0]]);
135  jacobian_vectors[1].set((2.0*vertices[v_v[2]]-vertices[v_v[0]]-
136  vertices[v_v[1]])*MSQ_SQRT_THREE_INV);
137  num_jacobian_vectors=2;
138  break;
139 
140  case QUADRILATERAL:
141  jacobian_vectors[0]=(vertices[v_v[1]]-vertices[v_v[0]]+sample_point[1]*
142  (vertices[v_v[2]]+vertices[v_v[0]]-vertices[v_v[3]]-
143  vertices[v_v[1]]));
144  jacobian_vectors[1]=(vertices[v_v[3]]-vertices[v_v[0]]+sample_point[0]*
145  (vertices[v_v[2]]+vertices[v_v[0]]-vertices[v_v[3]]-
146  vertices[v_v[1]]));
147  num_jacobian_vectors=2;
148  break;
149 
150  case TETRAHEDRON:
151  jacobian_vectors[0]=vertices[v_v[1]]-vertices[v_v[0]];
152  jacobian_vectors[1]=(2.0*vertices[v_v[2]]-vertices[v_v[0]]-
153  vertices[v_v[1]])*MSQ_SQRT_THREE_INV;
154  jacobian_vectors[2]=(3.0*vertices[v_v[3]]-vertices[v_v[2]]-
155  vertices[v_v[1]]-vertices[v_v[0]])*
156  MSQ_SQRT_TWO_INV*MSQ_SQRT_THREE_INV;
157  num_jacobian_vectors=3;
158  break;
159 
160  case HEXAHEDRON:
161 
162  jacobian_vectors[0]=vertices[v_v[1]]-vertices[v_v[0]]+
163  (sample_point[1]*(vertices[v_v[2]]+vertices[v_v[0]]-
164  vertices[v_v[3]]-vertices[v_v[1]]))+
165  (sample_point[2]*(vertices[v_v[5]]+vertices[v_v[0]]-
166  vertices[v_v[4]]-vertices[v_v[1]]))+
167  (sample_point[1]*sample_point[2]*(vertices[v_v[6]]+
168  vertices[v_v[4]]+
169  vertices[v_v[3]]+
170  vertices[v_v[1]]-
171  vertices[v_v[7]]-
172  vertices[v_v[5]]-
173  vertices[v_v[2]]-
174  vertices[v_v[0]])),
175  jacobian_vectors[1]=vertices[v_v[3]]-vertices[v_v[0]]+
176  (sample_point[0]*(vertices[v_v[2]]+vertices[v_v[0]]-
177  vertices[v_v[3]]-vertices[v_v[1]]))+
178  (sample_point[2]*(vertices[v_v[7]]+vertices[v_v[0]]-
179  vertices[v_v[4]]-vertices[v_v[3]]))+
180  (sample_point[0]*sample_point[2]*(vertices[v_v[6]]+
181  vertices[v_v[4]]+
182  vertices[v_v[3]]+
183  vertices[v_v[1]]-
184  vertices[v_v[7]]-
185  vertices[v_v[5]]-
186  vertices[v_v[2]]-
187  vertices[v_v[0]]));
188 
189  jacobian_vectors[2]=vertices[v_v[4]]-vertices[v_v[0]]+
190  (sample_point[0]*(vertices[v_v[5]]+vertices[v_v[0]]-
191  vertices[v_v[4]]-vertices[v_v[1]]))+
192  (sample_point[1]*(vertices[v_v[7]]+vertices[v_v[0]]-
193  vertices[v_v[4]]-vertices[v_v[3]]))+
194  (sample_point[0]*sample_point[1]*(vertices[v_v[6]]+
195  vertices[v_v[4]]+
196  vertices[v_v[3]]+
197  vertices[v_v[1]]-
198  vertices[v_v[7]]-
199  vertices[v_v[5]]-
200  vertices[v_v[2]]-
201  vertices[v_v[0]]));
202 
203  num_jacobian_vectors=3;
204  break;
205 
206  default:
207  MSQ_SETERR(err)(
208  "Compute_weighted_jacobian not yet defined for this entity.",
209  MsqError::NOT_IMPLEMENTED);
210  }
211 
212 }
213 
215 
222 void Mesquite::MsqMeshEntity::get_sample_points(QualityMetric::ElementEvaluationMode mode,
223  vector<Vector3D> &coords,
224  MsqError &err){
225  switch (mType)
226  {
227  case TRIANGLE:
228  switch (mode)
229  {
230  case (QualityMetric::ELEMENT_VERTICES):
231  coords.reserve(3);
232  coords.push_back(Vector3D(0.0, 0.0, 0.0));
233  coords.push_back(Vector3D(1.0, 0.0, 0.0));
234  coords.push_back(Vector3D(0.5, MSQ_SQRT_THREE_DIV_TWO, 0.0));
235  break;
236 
237  //The following need to be verified
238  case (QualityMetric::LINEAR_GAUSS_POINTS):
239  coords.reserve(1);
240  coords.push_back(Vector3D(0.5, MSQ_SQRT_THREE_DIV_TWO/2.0, 0.0));
241  break;
242 
243  case (QualityMetric::QUADRATIC_GAUSS_POINTS):
244  coords.reserve(3);
245  coords.push_back(Vector3D(0.5, 0.0, 0.0));
246  coords.push_back(Vector3D(0.75, MSQ_SQRT_THREE_DIV_TWO/2.0, 0.0));
247  coords.push_back(Vector3D(0.25, MSQ_SQRT_THREE_DIV_TWO/2.0, 0.0));
248  break;
249 
250  case (QualityMetric::CUBIC_GAUSS_POINTS):
251  coords.reserve(4);
252  coords.push_back(Vector3D(0.5, MSQ_SQRT_THREE_DIV_TWO/2.0, 0.0));
253  coords.push_back(Vector3D(0.2, 2.0* MSQ_SQRT_THREE_DIV_TWO/15.0,
254  0.0));
255  coords.push_back(Vector3D(0.8, 2.0* MSQ_SQRT_THREE_DIV_TWO/15.0,
256  0.0));
257  coords.push_back(Vector3D(0.5, 11.0* MSQ_SQRT_THREE_DIV_TWO/15.0,
258  0.0));
259  break;
260  default:
261  //return error saying sample points for mode not implem.
262  MSQ_SETERR(err)("Requested Sample Point Mode not implemented",
263  MsqError::NOT_IMPLEMENTED);
264  return;
265  }
266  break;
267 
268  case QUADRILATERAL:
269  switch (mode)
270  {
271  case (QualityMetric::ELEMENT_VERTICES):
272  coords.reserve(4);
273  coords.push_back(Vector3D(0.0, 0.0, 0.0));
274  coords.push_back(Vector3D(1.0, 0.0, 0.0));
275  coords.push_back(Vector3D(1.0, 1.0, 0.0));
276  coords.push_back(Vector3D(0.0, 1.0, 0.0));
277  break;
278  //THESE NEED TO BE VERIFIED
279  case (QualityMetric::LINEAR_GAUSS_POINTS):
280  coords.push_back(Vector3D(0.5, 0.5, 0.0));
281  break;
282  case (QualityMetric::QUADRATIC_GAUSS_POINTS):
283  case (QualityMetric::CUBIC_GAUSS_POINTS):
284  default:
285  //return error saying sample points for mode not implem.
286  MSQ_SETERR(err)("Requested Sample Point Mode not implemented",
287  MsqError::NOT_IMPLEMENTED);
288  return;
289  }
290  break;
291 
292  case TETRAHEDRON:
293  switch (mode)
294  {
295  case (QualityMetric::ELEMENT_VERTICES):
296  coords.reserve(4);
297  coords.push_back(Vector3D(0.0, 0.0, 0.0));
298  coords.push_back(Vector3D(1.0, 0.0, 0.0));
299  coords.push_back(Vector3D(0.5, MSQ_SQRT_THREE_DIV_TWO, 0.0));
300  coords.push_back(Vector3D(0.5, MSQ_SQRT_THREE_DIV_TWO/3.0,
301  MSQ_SQRT_TWO_DIV_SQRT_THREE));
302  break;
303  case (QualityMetric::LINEAR_GAUSS_POINTS):
304 
305  case (QualityMetric::QUADRATIC_GAUSS_POINTS):
306 
307  case (QualityMetric::CUBIC_GAUSS_POINTS):
308 
309  default:
310  //return error saying sample points for mode not implem.
311  MSQ_SETERR(err)("Requested Sample Point Mode not implemented",
312  MsqError::NOT_IMPLEMENTED);
313  return;
314  }
315  break;
316 
317  case HEXAHEDRON:
318  switch (mode)
319  {
320  case (QualityMetric::ELEMENT_VERTICES):
321  coords.reserve(8);
322  coords.push_back(Vector3D(0.0, 0.0, 0.0));
323  coords.push_back(Vector3D(1.0, 0.0, 0.0));
324  coords.push_back(Vector3D(1.0, 1.0, 0.0));
325  coords.push_back(Vector3D(0.0, 1.0, 0.0));
326  coords.push_back(Vector3D(0.0, 0.0, 1.0));
327  coords.push_back(Vector3D(1.0, 0.0, 1.0));
328  coords.push_back(Vector3D(1.0, 1.0, 1.0));
329  coords.push_back(Vector3D(0.0, 1.0, 1.0));
330  break;
331  case (QualityMetric::LINEAR_GAUSS_POINTS):
332 
333  case (QualityMetric::QUADRATIC_GAUSS_POINTS):
334 
335  case (QualityMetric::CUBIC_GAUSS_POINTS):
336 
337  default:
338  //return error saying sample points for mode not implem.
339  MSQ_SETERR(err)("Requested Sample Point Mode not implemented",
340  MsqError::NOT_IMPLEMENTED);
341  return;
342  }
343  break;
344 
345  default:
346  //return error saying sample points for mode not implem.
347  MSQ_SETERR(err)("Requested Sample Point Mode not implemented",
348  MsqError::NOT_IMPLEMENTED);
349  return;
350  }
351 }
352 
357 double MsqMeshEntity::compute_unsigned_area(PatchData &pd, MsqError &err) {
358  MsqVertex* verts=pd.get_vertex_array(err);MSQ_ERRZERO(err);
359  double tem=0.0;
360  switch (mType)
361  {
362 
363  case TRIANGLE:
364  tem = ((verts[vertexIndices[1]]-verts[vertexIndices[0]])*
365  (verts[vertexIndices[2]]-verts[vertexIndices[0]])).length()/2.0;
366  return tem;
367 
368  case QUADRILATERAL:
369  tem = ((verts[vertexIndices[1]]-verts[vertexIndices[0]])*
370  (verts[vertexIndices[3]]-verts[vertexIndices[0]])).length();
371  tem += ((verts[vertexIndices[3]]-verts[vertexIndices[2]])*
372  (verts[vertexIndices[1]]-verts[vertexIndices[2]])).length();
373  return (tem/2.0);
374 
375  default:
376  MSQ_SETERR(err)("Invalid type of element passed to compute unsigned area.",
377  MsqError::INVALID_ARG);
378  return 0;
379  }
380  return 0;
381 }
382 
387 double MsqMeshEntity::compute_unsigned_volume(PatchData &pd, MsqError &err) {
388  Vector3D sample_point(.5,.5,.5);
389  Vector3D jac_vecs[3];
390  short num_jacobian_vectors=-1;
391  double tem=0;
392  MsqVertex *verts = pd.get_vertex_array(err);MSQ_ERRZERO(err);
393  switch (mType)
394  {
395  case TETRAHEDRON:
396  tem = (verts[vertexIndices[3]]-verts[vertexIndices[0]])%
397  ((verts[vertexIndices[1]]-verts[vertexIndices[0]])*
398  (verts[vertexIndices[1]]-verts[vertexIndices[0]]))/6.0;
399  return fabs(tem);
400 
401  case HEXAHEDRON:
402  compute_weighted_jacobian(pd,sample_point,jac_vecs,
403  num_jacobian_vectors, err ); MSQ_ERRZERO(err);
404  return fabs(jac_vecs[2]%(jac_vecs[0]*jac_vecs[1]));
405 
406  default:
407  MSQ_SETERR(err)("Invalid type of element passed to compute unsigned volume.",
408  MsqError::INVALID_ARG);
409  return 0;
410  }
411  return 0;
412 }
413 
414 
419 double MsqMeshEntity::compute_signed_area(PatchData &pd, MsqError &err) {
420  MsqVertex* verts=pd.get_vertex_array(err);MSQ_ERRZERO(err);
421  double tem=0.0;
422  double tem2=0.0;
423  Vector3D surface_normal;
424  Vector3D cross_vec;
425  size_t element_index=pd.get_element_index(this);
426 
427  switch (mType)
428  {
429 
430  case TRIANGLE:
431  cross_vec=((verts[vertexIndices[1]]-verts[vertexIndices[0]])*
432  (verts[vertexIndices[2]]-verts[vertexIndices[0]]));
433  pd.get_domain_normal_at_element(element_index,surface_normal,err);
434  MSQ_ERRZERO(err);
435  tem = (cross_vec.length()/2.0);
436  //if normals do not point in same general direction, negate area
437  if(cross_vec%surface_normal<0){
438  tem *= -1;
439  }
440 
441  return tem;
442 
443  case QUADRILATERAL:
444  cross_vec=((verts[vertexIndices[1]]-verts[vertexIndices[0]])*
445  (verts[vertexIndices[3]]-verts[vertexIndices[0]]));
446  pd.get_domain_normal_at_element(element_index,surface_normal,err);
447  MSQ_ERRZERO(err);
448  tem = (cross_vec.length()/2.0);
449  //if normals do not point in same general direction, negate area
450  if(cross_vec%surface_normal<0){
451  tem *= -1;
452  }
453  cross_vec=((verts[vertexIndices[3]]-verts[vertexIndices[2]])*
454  (verts[vertexIndices[1]]-verts[vertexIndices[2]]));
455  tem2 = (cross_vec.length()/2.0);
456  //if normals do not point in same general direction, negate area
457  if(cross_vec%surface_normal<0){
458  tem2 *= -1;
459  //test to make sure surface normal existed
460  //if(surface_normal.length_squared()<.5){
461  //err.set_msg("compute_signed_area called without surface_normal available.");
462  //}
463  }
464  return (tem + tem2);
465 
466  default:
467  MSQ_SETERR(err)("Invalid type of element passed to compute unsigned area.",
468  MsqError::INVALID_ARG);
469  return 0;
470  };
471  return 0.0;
472 }
473 
478 double MsqMeshEntity::compute_signed_volume(PatchData &pd, MsqError &err) {
479  Vector3D sample_point(.5,.5,.5);
480  Vector3D jac_vecs[3];
481  short num_jacobian_vectors=-1;
482  double tem=0;
483  MsqVertex *verts = pd.get_vertex_array(err);MSQ_ERRZERO(err);
484  switch (mType)
485  {
486  case TETRAHEDRON:
487  tem = (verts[vertexIndices[3]]-verts[vertexIndices[0]])%
488  ((verts[vertexIndices[1]]-verts[vertexIndices[0]])*
489  (verts[vertexIndices[2]]-verts[vertexIndices[0]]))/6.0;
490  return tem;
491 
492  case HEXAHEDRON:
493  compute_weighted_jacobian(pd,sample_point,jac_vecs,
494  num_jacobian_vectors, err );
495  MSQ_ERRZERO(err);
496  return (jac_vecs[2]%(jac_vecs[0]*jac_vecs[1]));
497 
498  default:
499  MSQ_SETERR(err)("Invalid type of element passed to compute signed volume.",
500  MsqError::INVALID_ARG);
501  return 0;
502  };
503  return 0.0;
504 }
505 
506 
513 void Mesquite::MsqMeshEntity::get_connected_vertices(size_t vertex_index,
514  vector<size_t> &vert_indices,
515  MsqError &err)
516 {
517  //i iterates through elem's vertices
518  int i=0;
519  //index is set to the index in the vertexIndices corresponding
520  //to vertex_index
521  int index=-1;
522 
523  switch (mType)
524  {
525  case TRIANGLE:
526  while(i<3)
527  {
528  if(vertexIndices[i]==vertex_index)
529  {
530  index=i;
531  break;
532  }
533  ++i;
534  }
535  if(index>=0)
536  {
537  vert_indices.push_back(vertexIndices[(index+1)%3]);
538  vert_indices.push_back(vertexIndices[(index+2)%3]);
539  }
540 
541  break;
542 
543  case QUADRILATERAL:
544  while(i<4)
545  {
546  if(vertexIndices[i]==vertex_index)
547  {
548  index=i;
549  break;
550  }
551  ++i;
552  }
553  if(index>=0)
554  {
555  vert_indices.push_back(vertexIndices[(index+1)%4]);
556  vert_indices.push_back(vertexIndices[(index+3)%4]);
557  }
558 
559  break;
560 
561  case TETRAHEDRON:
562  while(i<4)
563  {
564  if(vertexIndices[i]==vertex_index)
565  {
566  index=i;
567  break;
568  }
569  ++i;
570  }
571  if(index>=0)
572  {
573  vert_indices.push_back(vertexIndices[(index+1)%4]);
574  vert_indices.push_back(vertexIndices[(index+2)%4]);
575  vert_indices.push_back(vertexIndices[(index+3)%4]);
576  }
577 
578  break;
579 
580  case HEXAHEDRON:
581  while(i<8)
582  {
583  if(vertexIndices[i]==vertex_index)
584  {
585  index=i;
586  break;
587  }
588  ++i;
589  }
590 
591  if(index>=0)
592  {
593  if (index<4)
594  {
595  vert_indices.push_back(vertexIndices[(index+1)%4]);
596  vert_indices.push_back(vertexIndices[(index+3)%4]);
597  vert_indices.push_back(vertexIndices[(index)+4]);
598  }
599  else
600  {
601  vert_indices.push_back(vertexIndices[(index+3)%4+4]);
602  vert_indices.push_back(vertexIndices[(index+1)%4+4]);
603  vert_indices.push_back(vertexIndices[(index)-4]);
604  }
605  }
606 
607  break;
608 
609  default:
610  MSQ_SETERR(err)("Element type not available", MsqError::INVALID_ARG);
611  break;
612  }
613 }
614 
618 /*
619 void MsqMeshEntity::compute_corner_normal(size_t corner,
620  Vector3D &normal,
621  PatchData &pd,
622  MsqError &err)
623 {
624  if ( get_element_type()==TRIANGLE || get_element_type()==QUADRILATERAL )
625  {
626  // There are two cases where we cannot get a normal from the
627  // geometry that are not errors:
628  // 1) There is no domain set
629  // 2) The vertex is at a degenerate point on the geometry (e.g.
630  // tip of a cone.)
631 
632  bool have_normal = false;
633 
634  // Get normal from domain
635  if (pd.domain_set())
636  {
637  size_t index = pd.get_element_index(this);
638  pd.get_domain_normal_at_corner( index, corner, normal, err );
639  MSQ_ERRRTN(err);
640 
641  double length = normal.length();
642  if (length > DBL_EPSILON)
643  {
644  have_normal = true;
645  normal /= length;
646  }
647  }
648 
649  // If failed to get normal from domain, calculate
650  // from adjacent vertices.
651  if ( !have_normal )
652  {
653  const int num_corner = this->vertex_count();
654  const int prev_corner = (corner + num_corner - 1) % num_corner;
655  const int next_corner = (corner + 1 ) % num_corner;
656  const size_t this_idx = vertexIndices[corner];
657  const size_t prev_idx = vertexIndices[prev_corner];
658  const size_t next_idx = vertexIndices[next_corner];
659  normal = (pd.vertex_by_index(next_idx) - pd.vertex_by_index(this_idx))
660  * (pd.vertex_by_index(prev_idx) - pd.vertex_by_index(this_idx));
661  normal.normalize();
662  }
663  }
664  else
665  MSQ_SETERR(err)("Should only be used for faces (tri, quads, ...).",
666  MsqError::INVALID_ARG);
667 }
668 */
669 
670 void MsqMeshEntity::compute_corner_normals( Vector3D normals[],
671  PatchData &pd,
672  MsqError &err)
673 {
674  EntityTopology type = get_element_type();
675  if (type != TRIANGLE && type != QUADRILATERAL && type != POLYGON)
676  {
677  MSQ_SETERR(err)("Should only be used for faces (tri, quads, ...).",
678  MsqError::INVALID_ARG);
679  return;
680  }
681 
682 
683  // There are two cases where we cannot get a normal from the
684  // geometry that are not errors:
685  // 1) There is no domain set
686  // 2) The vertex is at a degenerate point on the geometry (e.g.
687  // tip of a cone.)
688 
689  // Get normal from domain
690  if (pd.domain_set())
691  {
692  size_t index = pd.get_element_index(this);
693  pd.get_domain_normals_at_corners( index, normals, err );
694  MSQ_ERRRTN(err);
695  }
696 
697  // Check if normals are valid (none are valid if !pd.domain_set())
698  const unsigned count = vertex_count();
699  for (unsigned i = 0; i < count; ++i)
700  {
701  // If got valid normal from domain,
702  // make it a unit vector and continue.
703  if (pd.domain_set())
704  {
705  double length = normals[i].length();
706  if (length > DBL_EPSILON)
707  {
708  normals[i] /= length;
709  continue;
710  }
711  }
712 
713  const size_t prev_idx = vertexIndices[(i + count - 1) % count];
714  const size_t this_idx = vertexIndices[i];
715  const size_t next_idx = vertexIndices[(i + 1) % count];
716 
717  // Calculate normal using edges adjacent to corner
718  normals[i] = (pd.vertex_by_index(next_idx) - pd.vertex_by_index(this_idx))
719  * (pd.vertex_by_index(prev_idx) - pd.vertex_by_index(this_idx));
720  normals[i].normalize();
721  }
722 }
723 
731 void MsqMeshEntity::compute_corner_matrices(PatchData &pd, Matrix3D A[], int num_m3d, MsqError &err )
732 {
733  MsqVertex* vertices = pd.get_vertex_array(err); MSQ_ERRRTN(err);
734  const size_t* v_i = &vertexIndices[0];
735 
736  // If 2D element, we will get the surface normal
737  Vector3D normals[4], vec1, vec2, vec3, vec4;
738 
739 
740  switch(get_element_type()){
741 
742  case TRIANGLE:
743  if (num_m3d != 3) {
744  MSQ_SETERR(err)("num_m3d incompatible with element type.", MsqError::INVALID_ARG);
745  return;
746  }
747 
748  compute_corner_normals( normals, pd, err ); MSQ_ERRRTN(err);
749 
750  vec1 = vertices[v_i[1]]-vertices[v_i[0]];
751  vec2 = vertices[v_i[2]]-vertices[v_i[0]];
752  vec3 = vertices[v_i[2]]-vertices[v_i[1]];
753 
754  A[0].set_column(0, vec1);
755  A[0].set_column(1, vec2);
756  A[0].set_column(2, normals[0]*MSQ_3RT_2_OVER_6RT_3);
757 
758  A[1].set_column(0, vec3);
759  A[1].set_column(1, -vec1);
760  A[1].set_column(2, normals[1]*MSQ_3RT_2_OVER_6RT_3);
761 
762  A[2].set_column(0, -vec2);
763  A[2].set_column(1, -vec3);
764  A[2].set_column(2, normals[2]*MSQ_3RT_2_OVER_6RT_3);
765 
766  break;
767 
768  case QUADRILATERAL:
769  if (num_m3d != 4) {
770  MSQ_SETERR(err)("num_m3d incompatible with element type.", MsqError::INVALID_ARG);
771  return;
772  }
773  vec1 = vertices[v_i[1]]-vertices[v_i[0]];
774  vec2 = vertices[v_i[3]]-vertices[v_i[0]];
775  vec3 = vertices[v_i[2]]-vertices[v_i[1]];
776  vec4 = vertices[v_i[3]]-vertices[v_i[2]];
777 
778  compute_corner_normals( normals, pd, err ); MSQ_ERRRTN(err);
779 
780  A[0].set_column(0, vec1);
781  A[0].set_column(1, vec2);
782  A[0].set_column(2, normals[0]);
783 
784  A[1].set_column(0, vec3);
785  A[1].set_column(1, -vec1);
786  A[1].set_column(2, normals[1]);
787 
788  A[2].set_column(0, vec4);
789  A[2].set_column(1, -vec3);
790  A[2].set_column(2, normals[2]);
791 
792  A[3].set_column(0, -vec2);
793  A[3].set_column(1, -vec4);
794  A[3].set_column(2, normals[3]);
795 
796  break;
797 
798  case TETRAHEDRON:
799  if (num_m3d != 4) {
800  MSQ_SETERR(err)("num_m3d incompatible with element type.", MsqError::INVALID_ARG);
801  return;
802  }
803  A[0].set_column(0, vertices[v_i[1]]-vertices[v_i[0]]);
804  A[0].set_column(1, vertices[v_i[2]]-vertices[v_i[0]]);
805  A[0].set_column(2, vertices[v_i[3]]-vertices[v_i[0]]);
806 
807  A[1].set_column(0, vertices[v_i[0]]-vertices[v_i[1]]);
808  A[1].set_column(1, vertices[v_i[3]]-vertices[v_i[1]]);
809  A[1].set_column(2, vertices[v_i[2]]-vertices[v_i[1]]);
810 
811  A[2].set_column(0, vertices[v_i[3]]-vertices[v_i[2]]);
812  A[2].set_column(1, vertices[v_i[0]]-vertices[v_i[2]]);
813  A[2].set_column(2, vertices[v_i[1]]-vertices[v_i[2]]);
814 
815  A[3].set_column(0, vertices[v_i[2]]-vertices[v_i[3]]);
816  A[3].set_column(1, vertices[v_i[1]]-vertices[v_i[3]]);
817  A[3].set_column(2, vertices[v_i[0]]-vertices[v_i[3]]);
818 
819  break;
820 
821  /*
822  case PYRAMID:
823  //We compute the pyramid's "condition number" by averaging
824  //the 4 tet's condition numbers, where the tets are created
825  //by removing one of the four base vertices from the pyramid.
826  //transform to origina v_i[0]
827  temp_vec[3]=vertices[v_i[1]]-vertices[v_i[0]];
828  temp_vec[4]=vertices[v_i[3]]-vertices[v_i[0]];
829  temp_vec[5]=vertices[v_i[4]]-vertices[v_i[0]];
830  //find AW_inverse
831  temp_vec[0]=temp_vec[3];
832  temp_vec[1]=temp_vec[4]-temp_vec[3];
833  temp_vec[2]=MSQ_SQRT_TWO*(temp_vec[5]-(temp_vec[4]/2.0));
834  return_flag=condition_number_3d(temp_vec,pd,met_vals[0],err);
835  if(!return_flag)
836  return return_flag;
837  //transform to origina v_i[1]
838  temp_vec[3]=vertices[v_i[2]]-vertices[v_i[1]];
839  temp_vec[4]=vertices[v_i[3]]-vertices[v_i[1]];
840  temp_vec[5]=vertices[v_i[4]]-vertices[v_i[1]];
841  //find AW_inverse
842  temp_vec[0]=temp_vec[3]-temp_vec[4];
843  temp_vec[1]=temp_vec[3];
844  temp_vec[2]=MSQ_SQRT_TWO*(temp_vec[5]-(temp_vec[4]/2.0));
845  return_flag=condition_number_3d(temp_vec,pd,met_vals[1],err);
846  if(!return_flag)
847  return return_flag;
848  //transform to origina v_i[1]
849  temp_vec[3]=vertices[v_i[3]]-vertices[v_i[2]];
850  temp_vec[4]=vertices[v_i[0]]-vertices[v_i[2]];
851  temp_vec[5]=vertices[v_i[4]]-vertices[v_i[2]];
852  //find AW_inverse
853  temp_vec[0]=-temp_vec[3];
854  temp_vec[1]=temp_vec[3]-temp_vec[4];
855  temp_vec[2]=MSQ_SQRT_TWO*(temp_vec[5]-(temp_vec[4]/2.0));
856  return_flag=condition_number_3d(temp_vec,pd,met_vals[2],err);
857  if(!return_flag)
858  return return_flag;
859  //transform to origina v_i[1]
860  temp_vec[3]=vertices[v_i[0]]-vertices[v_i[3]];
861  temp_vec[4]=vertices[v_i[1]]-vertices[v_i[3]];
862  temp_vec[5]=vertices[v_i[4]]-vertices[v_i[3]];
863  //find AW_inverse
864  temp_vec[0]=temp_vec[4]-temp_vec[3];
865  temp_vec[1]=-temp_vec[3];
866  temp_vec[2]=MSQ_SQRT_TWO*(temp_vec[5]-(temp_vec[4]/2.0));
867  return_flag=condition_number_3d(temp_vec,pd,met_vals[3],err);
868  fval=average_metrics(met_vals, 4, err);
869  if(!return_flag)
870  return return_flag;
871  break;
872  */
873 
874  case HEXAHEDRON:
875  if (num_m3d != 8) {
876  MSQ_SETERR(err)("num_m3d incompatible with element type.", MsqError::INVALID_ARG);
877  return;
878  }
879  A[0].set_column(0, vertices[v_i[1]]-vertices[v_i[0]]);
880  A[0].set_column(1, vertices[v_i[3]]-vertices[v_i[0]]);
881  A[0].set_column(2, vertices[v_i[4]]-vertices[v_i[0]]);
882 
883  A[1].set_column(0, vertices[v_i[2]]-vertices[v_i[1]]);
884  A[1].set_column(1, vertices[v_i[0]]-vertices[v_i[1]]);
885  A[1].set_column(2, vertices[v_i[5]]-vertices[v_i[1]]);
886 
887  A[2].set_column(0, vertices[v_i[3]]-vertices[v_i[2]]);
888  A[2].set_column(1, vertices[v_i[1]]-vertices[v_i[2]]);
889  A[2].set_column(2, vertices[v_i[6]]-vertices[v_i[2]]);
890 
891  A[3].set_column(0, vertices[v_i[0]]-vertices[v_i[3]]);
892  A[3].set_column(1, vertices[v_i[2]]-vertices[v_i[3]]);
893  A[3].set_column(2, vertices[v_i[7]]-vertices[v_i[3]]);
894 
895  A[4].set_column(0, vertices[v_i[7]]-vertices[v_i[4]]);
896  A[4].set_column(1, vertices[v_i[5]]-vertices[v_i[4]]);
897  A[4].set_column(2, vertices[v_i[0]]-vertices[v_i[4]]);
898 
899  A[5].set_column(0, vertices[v_i[4]]-vertices[v_i[5]]);
900  A[5].set_column(1, vertices[v_i[6]]-vertices[v_i[5]]);
901  A[5].set_column(2, vertices[v_i[1]]-vertices[v_i[5]]);
902 
903  A[6].set_column(0, vertices[v_i[5]]-vertices[v_i[6]]);
904  A[6].set_column(1, vertices[v_i[7]]-vertices[v_i[6]]);
905  A[6].set_column(2, vertices[v_i[2]]-vertices[v_i[6]]);
906 
907  A[7].set_column(0, vertices[v_i[6]]-vertices[v_i[7]]);
908  A[7].set_column(1, vertices[v_i[4]]-vertices[v_i[7]]);
909  A[7].set_column(2, vertices[v_i[3]]-vertices[v_i[7]]);
910 
911  break;
912 
913  default:
914  MSQ_SETERR(err)("element type not implemented.", MsqError::NOT_IMPLEMENTED);
915  return;
916  }// end switch over element type
917 }
918 
919 ostream& operator<<( ostream& stream, const MsqMeshEntity& entity )
920 {
921  size_t num_vert = entity.vertex_count();
922  stream << "MsqMeshEntity " << &entity << " with vertices ";
923  for (size_t i = 0; i < num_vert; ++i)
924  stream << entity.vertexIndices[i] << " ";
925  stream << endl;
926  return stream;
927 }
928 
929 } // namespace Mesquite
MSQ_ERRZERO
#define MSQ_ERRZERO(err)
Return zero/NULL on error.
Definition: includeLinks/MsqError.hpp:72
Mesquite::MsqMeshEntity::get_node_indices
void get_node_indices(msq_std::vector< msq_stdc::size_t > &vertex_list) const
gets the vertices of the mesh entity
Mesquite::EntityTopology
EntityTopology
Definition: Mesquite.hpp:92
Mesquite::MsqMeshEntity::compute_signed_area
double compute_signed_area(PatchData &pd, MsqError &err)
Computes the signed area of the element.
Definition: Mesh/MsqMeshEntity.cpp:419
Mesquite::MsqError::NOT_IMPLEMENTED
requested functionality is not (yet) implemented
Definition: includeLinks/MsqError.hpp:128
Mesquite::MsqMeshEntity::get_vertex_indices
void get_vertex_indices(msq_std::vector< msq_stdc::size_t > &vertex_list) const
gets the vertices of the mesh entity
Mesquite::MSQ_SQRT_TWO_DIV_SQRT_THREE
static const double MSQ_SQRT_TWO_DIV_SQRT_THREE
Definition: Mesquite.hpp:127
Mesquite::MsqMeshEntity::compute_unsigned_volume
double compute_unsigned_volume(PatchData &pd, MsqError &err)
Computes the volume of the element.
Definition: Mesh/MsqMeshEntity.cpp:387
Mesquite::MSQ_SQRT_THREE_DIV_TWO
static const double MSQ_SQRT_THREE_DIV_TWO
Definition: Mesquite.hpp:124
Mesquite::MsqMeshEntity::vertexIndices
size_t * vertexIndices
Pointer to connectivity array.
Definition: includeLinks/MsqMeshEntity.hpp:173
Mesquite::length
double length(Vector3D *const v, int n)
Definition: includeLinks/Vector3D.hpp:400
Mesquite::MsqError::INVALID_ARG
invalid function argument passed
Definition: includeLinks/MsqError.hpp:119
A
rational * A
Definition: vinci_lass.c:67
nvc::Vector3D
NVec< 3, double > Vector3D
Definition: Rocon/include/NVec.h:219
Mesquite::MsqMeshEntity::vertex_count
msq_stdc::size_t vertex_count() const
Returns the number of vertices in this element, based on its element type.
Definition: includeLinks/MsqMeshEntity.hpp:185
MSQ_SETERR
#define MSQ_SETERR(err)
Macro to set error - use err.clear() to clear.
Definition: includeLinks/MsqError.hpp:83
i
blockLoc i
Definition: read.cpp:79
Mesquite::MsqMeshEntity::append_vertex_indices
void append_vertex_indices(msq_std::vector< msq_stdc::size_t > &vertex_list) const
Mesquite::MsqMeshEntity::compute_weighted_jacobian
void compute_weighted_jacobian(PatchData &pd, Vector3D &sample_point, Vector3D jacobian_vectors[], short &num_jacobian_vectors, MsqError &err)
Definition: Mesh/MsqMeshEntity.cpp:118
Mesquite::MsqMeshEntity::get_centroid
void get_centroid(Vector3D &centroid, const PatchData &pd, MsqError &err) const
Returns the centroid of the element.
Definition: Mesh/MsqMeshEntity.cpp:106
Mesquite::MsqMeshEntity::compute_unsigned_area
double compute_unsigned_area(PatchData &pd, MsqError &err)
Computes the area of the element.
Definition: Mesh/MsqMeshEntity.cpp:357
Mesquite::MsqMeshEntity::append_node_indices
void append_node_indices(msq_std::vector< msq_stdc::size_t > &vertex_list) const
Mesquite::MSQ_SQRT_TWO_INV
static const double MSQ_SQRT_TWO_INV
Definition: Mesquite.hpp:126
Mesquite::MsqMeshEntity::get_connected_vertices
void get_connected_vertices(msq_stdc::size_t vertex_index, msq_std::vector< msq_stdc::size_t > &vert_indices, MsqError &err)
Fills a vector&lt;size_t&gt; with vertices connected to the given vertex through the edges of this MsqMeshE...
Mesquite::MsqMeshEntity::get_element_type
EntityTopology get_element_type() const
Returns element type.
Definition: includeLinks/MsqMeshEntity.hpp:78
Mesquite::MsqMeshEntity::compute_signed_volume
double compute_signed_volume(PatchData &pd, MsqError &err)
Computes the signed volume of the element.
Definition: Mesh/MsqMeshEntity.cpp:478
Mesquite::MSQ_SQRT_THREE_INV
static const double MSQ_SQRT_THREE_INV
Definition: Mesquite.hpp:125
Mesquite::operator<<
msq_stdio::ostream & operator<<(msq_stdio::ostream &s, const Matrix3D &A)
Definition: includeLinks/Matrix3D.hpp:250
Mesquite::MsqMeshEntity::get_sample_points
void get_sample_points(QualityMetric::ElementEvaluationMode mode, msq_std::vector< Vector3D > &coords, MsqError &err)
Returns a list of sample points given an evaluationmode.
Mesquite::MsqMeshEntity::compute_corner_matrices
void compute_corner_matrices(PatchData &pd, Matrix3D A[], int num_m3d, MsqError &err)
Compute matrices which column are the vectors issued from a corner.
Definition: Mesh/MsqMeshEntity.cpp:731
MSQ_ERRRTN
#define MSQ_ERRRTN(err)
If passed error is true, return from a void function.
Definition: includeLinks/MsqError.hpp:68
Mesquite::MsqMeshEntity::compute_corner_normals
void compute_corner_normals(Vector3D normals[], PatchData &pd, MsqError &err)
Uses a MeshDomain call-back function to compute the normal at the corner.
Definition: Mesh/MsqMeshEntity.cpp:670
Mesquite::MSQ_3RT_2_OVER_6RT_3
static const double MSQ_3RT_2_OVER_6RT_3
Definition: Mesquite.hpp:130
Mesquite::MsqMeshEntity::node_count
msq_stdc::size_t node_count() const
Return number of nodes in element (number of corner vertices + number of higher-order nodes)...
Definition: includeLinks/MsqMeshEntity.hpp:86