Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
inks/ConditionNumberQualityMetric.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  ***************************************************************** */
34 #include <vector>
35 #include "ConditionNumberQualityMetric.hpp"
36 #include <math.h>
37 #include "Vector3D.hpp"
38 #include "ShapeQualityMetric.hpp"
39 #include "QualityMetric.hpp"
40 
41 using namespace Mesquite;
42 
43 ConditionNumberQualityMetric::ConditionNumberQualityMetric()
44 {
45  MsqError err;
46  set_metric_type(ELEMENT_BASED);
47  set_element_evaluation_mode(ELEMENT_VERTICES, err); MSQ_CHKERR(err);
48  avgMethod=QualityMetric::LINEAR;
49  feasible=1;
50  set_name("Condition Number");
51 }
52 
53 
54 bool ConditionNumberQualityMetric::evaluate_element(PatchData &pd,
55  MsqMeshEntity *element,
56  double &fval,
57  MsqError &err)
58 {
59  bool return_flag;
60  double met_vals[MSQ_MAX_NUM_VERT_PER_ENT];
61  fval=MSQ_MAX_CAP;
62  const size_t* v_i = element->get_vertex_index_array();
63  //only 3 temp_vec will be sent to cond-num calculator, but the
64  //additional vector3Ds may be needed during the calculations
65  Vector3D temp_vec[6];
66  MsqVertex *vertices=pd.get_vertex_array(err);
67  switch(element->get_element_type()){
68  case TRIANGLE:
69  temp_vec[0]=vertices[v_i[1]]-vertices[v_i[0]];
70  temp_vec[2]=vertices[v_i[2]]-vertices[v_i[0]];
71  //make relative to equilateral
72  temp_vec[1]=((2*temp_vec[2])-temp_vec[0])*MSQ_SQRT_THREE_INV;
73  return_flag=condition_number_2d(temp_vec,v_i[0],pd,fval,err); MSQ_ERRZERO(err);
74  return return_flag;
75  case QUADRILATERAL:
76  temp_vec[0]=vertices[v_i[1]]-vertices[v_i[0]];
77  temp_vec[1]=vertices[v_i[3]]-vertices[v_i[0]];
78  return_flag=condition_number_2d(temp_vec,v_i[0],pd,met_vals[0],err); MSQ_ERRZERO(err);
79  if(!return_flag)
80  return return_flag;
81  temp_vec[0]=vertices[v_i[2]]-vertices[v_i[1]];
82  temp_vec[1]=vertices[v_i[0]]-vertices[v_i[1]];
83  return_flag=condition_number_2d(temp_vec,v_i[1],pd,met_vals[1],err); MSQ_ERRZERO(err);
84  if(!return_flag)
85  return return_flag;
86  temp_vec[0]=vertices[v_i[3]]-vertices[v_i[2]];
87  temp_vec[1]=vertices[v_i[1]]-vertices[v_i[2]];
88  return_flag=condition_number_2d(temp_vec,v_i[2],pd,met_vals[2],err); MSQ_ERRZERO(err);
89  if(!return_flag)
90  return return_flag;
91  temp_vec[0]=vertices[v_i[0]]-vertices[v_i[3]];
92  temp_vec[1]=vertices[v_i[2]]-vertices[v_i[3]];
93  return_flag=condition_number_2d(temp_vec,v_i[3],pd,met_vals[3],err); MSQ_ERRZERO(err);
94  fval = average_metrics(met_vals, 4, err);
95  return return_flag;
96  case TETRAHEDRON:
97  temp_vec[0]=vertices[v_i[1]]-vertices[v_i[0]];
98  temp_vec[3]=vertices[v_i[2]]-vertices[v_i[0]];
99  temp_vec[4]=vertices[v_i[3]]-vertices[v_i[0]];
100  //transform to equilateral tet
101  temp_vec[1]=((2*temp_vec[3])-temp_vec[0])/MSQ_SQRT_THREE;
102  temp_vec[2]=((3*temp_vec[4])-temp_vec[0]-temp_vec[3])/
103  (MSQ_SQRT_THREE*MSQ_SQRT_TWO);
104  return_flag=condition_number_3d(temp_vec,pd,fval,err); MSQ_ERRZERO(err);
105  return return_flag;
106  /*
107  case PYRAMID:
108  //We compute the pyramid's "condition number" by averaging
109  //the 4 tet's condition numbers, where the tets are created
110  //by removing one of the four base vertices from the pyramid.
111  //transform to origina v_i[0]
112  temp_vec[3]=vertices[v_i[1]]-vertices[v_i[0]];
113  temp_vec[4]=vertices[v_i[3]]-vertices[v_i[0]];
114  temp_vec[5]=vertices[v_i[4]]-vertices[v_i[0]];
115  //find AW_inverse
116  temp_vec[0]=temp_vec[3];
117  temp_vec[1]=temp_vec[4]-temp_vec[3];
118  temp_vec[2]=MSQ_SQRT_TWO*(temp_vec[5]-(temp_vec[4]/2.0));
119  return_flag=condition_number_3d(temp_vec,pd,met_vals[0],err);
120  if(!return_flag)
121  return return_flag;
122  //transform to origina v_i[1]
123  temp_vec[3]=vertices[v_i[2]]-vertices[v_i[1]];
124  temp_vec[4]=vertices[v_i[3]]-vertices[v_i[1]];
125  temp_vec[5]=vertices[v_i[4]]-vertices[v_i[1]];
126  //find AW_inverse
127  temp_vec[0]=temp_vec[3]-temp_vec[4];
128  temp_vec[1]=temp_vec[3];
129  temp_vec[2]=MSQ_SQRT_TWO*(temp_vec[5]-(temp_vec[4]/2.0));
130  return_flag=condition_number_3d(temp_vec,pd,met_vals[1],err);
131  if(!return_flag)
132  return return_flag;
133  //transform to origina v_i[1]
134  temp_vec[3]=vertices[v_i[3]]-vertices[v_i[2]];
135  temp_vec[4]=vertices[v_i[0]]-vertices[v_i[2]];
136  temp_vec[5]=vertices[v_i[4]]-vertices[v_i[2]];
137  //find AW_inverse
138  temp_vec[0]=-temp_vec[3];
139  temp_vec[1]=temp_vec[3]-temp_vec[4];
140  temp_vec[2]=MSQ_SQRT_TWO*(temp_vec[5]-(temp_vec[4]/2.0));
141  return_flag=condition_number_3d(temp_vec,pd,met_vals[2],err);
142  if(!return_flag)
143  return return_flag;
144  //transform to origina v_i[1]
145  temp_vec[3]=vertices[v_i[0]]-vertices[v_i[3]];
146  temp_vec[4]=vertices[v_i[1]]-vertices[v_i[3]];
147  temp_vec[5]=vertices[v_i[4]]-vertices[v_i[3]];
148  //find AW_inverse
149  temp_vec[0]=temp_vec[4]-temp_vec[3];
150  temp_vec[1]=-temp_vec[3];
151  temp_vec[2]=MSQ_SQRT_TWO*(temp_vec[5]-(temp_vec[4]/2.0));
152  return_flag=condition_number_3d(temp_vec,pd,met_vals[3],err);
153  fval=average_metrics(met_vals, 4, err);
154  if(!return_flag)
155  return return_flag;
156  break;
157  */
158  case HEXAHEDRON:
159  //transform to v_i[0]
160  temp_vec[0]=vertices[v_i[1]]-vertices[v_i[0]];
161  temp_vec[1]=vertices[v_i[3]]-vertices[v_i[0]];
162  temp_vec[2]=vertices[v_i[4]]-vertices[v_i[0]];
163  return_flag=condition_number_3d(temp_vec,pd,met_vals[0],err); MSQ_ERRZERO(err);
164  if(!return_flag)
165  return return_flag;
166  temp_vec[0]=vertices[v_i[2]]-vertices[v_i[1]];
167  temp_vec[1]=vertices[v_i[0]]-vertices[v_i[1]];
168  temp_vec[2]=vertices[v_i[5]]-vertices[v_i[1]];
169  return_flag=condition_number_3d(temp_vec,pd,met_vals[1],err); MSQ_ERRZERO(err);
170  if(!return_flag)
171  return return_flag;
172  temp_vec[0]=vertices[v_i[3]]-vertices[v_i[2]];
173  temp_vec[1]=vertices[v_i[1]]-vertices[v_i[2]];
174  temp_vec[2]=vertices[v_i[6]]-vertices[v_i[2]];
175  return_flag=condition_number_3d(temp_vec,pd,met_vals[2],err); MSQ_ERRZERO(err);
176  if(!return_flag)
177  return return_flag;
178  temp_vec[0]=vertices[v_i[0]]-vertices[v_i[3]];
179  temp_vec[1]=vertices[v_i[2]]-vertices[v_i[3]];
180  temp_vec[2]=vertices[v_i[7]]-vertices[v_i[3]];
181  return_flag=condition_number_3d(temp_vec,pd,met_vals[3],err); MSQ_ERRZERO(err);
182  if(!return_flag)
183  return return_flag;
184  temp_vec[0]=vertices[v_i[7]]-vertices[v_i[4]];
185  temp_vec[1]=vertices[v_i[5]]-vertices[v_i[4]];
186  temp_vec[2]=vertices[v_i[0]]-vertices[v_i[4]];
187  return_flag=condition_number_3d(temp_vec,pd,met_vals[4],err); MSQ_ERRZERO(err);
188  if(!return_flag)
189  return return_flag;
190  temp_vec[0]=vertices[v_i[4]]-vertices[v_i[5]];
191  temp_vec[1]=vertices[v_i[6]]-vertices[v_i[5]];
192  temp_vec[2]=vertices[v_i[1]]-vertices[v_i[5]];
193  return_flag=condition_number_3d(temp_vec,pd,met_vals[5],err); MSQ_ERRZERO(err);
194  if(!return_flag)
195  return return_flag;
196  temp_vec[0]=vertices[v_i[5]]-vertices[v_i[6]];
197  temp_vec[1]=vertices[v_i[7]]-vertices[v_i[6]];
198  temp_vec[2]=vertices[v_i[2]]-vertices[v_i[6]];
199  return_flag=condition_number_3d(temp_vec,pd,met_vals[6],err); MSQ_ERRZERO(err);
200  if(!return_flag)
201  return return_flag;
202  temp_vec[0]=vertices[v_i[6]]-vertices[v_i[7]];
203  temp_vec[1]=vertices[v_i[4]]-vertices[v_i[7]];
204  temp_vec[2]=vertices[v_i[3]]-vertices[v_i[7]];
205  return_flag=condition_number_3d(temp_vec,pd,met_vals[7],err); MSQ_ERRZERO(err);
206  fval=average_metrics(met_vals, 8, err); MSQ_ERRZERO(err);
207  return return_flag;
208  default:
209  fval=MSQ_MAX_CAP;
210  }// end switch over element type
211  return false;
212 }
213 
214 
Mesquite::ShapeQualityMetric::condition_number_3d
bool condition_number_3d(Vector3D temp_vec[], PatchData &pd, double &fval, MsqError &err)
Given the 3-d jacobian matrix, compute the condition number, fval.
Definition: includeLinks/ShapeQualityMetric.hpp:186
MSQ_ERRZERO
#define MSQ_ERRZERO(err)
Return zero/NULL on error.
Definition: includeLinks/MsqError.hpp:72
Mesquite::MSQ_MAX_CAP
const double MSQ_MAX_CAP
Definition: Mesquite.hpp:173
Mesquite::MsqError
Used to hold the error state and return it to the application.
Definition: includeLinks/MsqError.hpp:106
Mesquite::ConditionNumberQualityMetric::evaluate_element
bool evaluate_element(PatchData &pd, MsqMeshEntity *element, double &fval, MsqError &err)
evaluate using mesquite objects
Definition: QualityMetric/Shape/ConditionNumberQualityMetric.cpp:54
Mesquite::MsqMeshEntity
MsqMeshEntity is the Mesquite object that stores information about the elements in the mesh...
Definition: includeLinks/MsqMeshEntity.hpp:69
Mesquite::Vector3D
Vector3D is the object that effeciently stores information about about three-deminsional vectors...
Definition: includeLinks/Vector3D.hpp:64
Mesquite::MSQ_MAX_NUM_VERT_PER_ENT
const int MSQ_MAX_NUM_VERT_PER_ENT
Definition: Mesquite.hpp:120
Mesquite::MSQ_SQRT_TWO
static const double MSQ_SQRT_TWO
Definition: Mesquite.hpp:122
Mesquite::QualityMetric::set_element_evaluation_mode
void set_element_evaluation_mode(ElementEvaluationMode mode, MsqError &err)
Sets the evaluation mode for the ELEMENT_BASED metrics.
Definition: includeLinks/QualityMetric.hpp:402
MSQ_CHKERR
#define MSQ_CHKERR(err)
Mesquite&#39;s Error Checking macro.
Definition: includeLinks/MsqError.hpp:62
Mesquite::ShapeQualityMetric::condition_number_2d
bool condition_number_2d(Vector3D temp_vec[], size_t v_ind, PatchData &pd, double &fval, MsqError &err)
Given the 2-d jacobian matrix, compute the condition number, fval.
Definition: includeLinks/ShapeQualityMetric.hpp:109
Mesquite::PatchData::get_vertex_array
const MsqVertex * get_vertex_array(MsqError &err) const
Returns a pointer to the start of the vertex array.
Definition: includeLinks/PatchData.hpp:513
Mesquite::MsqMeshEntity::get_element_type
EntityTopology get_element_type() const
Returns element type.
Definition: includeLinks/MsqMeshEntity.hpp:78
Mesquite::MsqMeshEntity::get_vertex_index_array
const msq_stdc::size_t * get_vertex_index_array() const
Very efficient retrieval of vertices indexes (corresponding to the PatchData vertex array)...
Definition: includeLinks/MsqMeshEntity.hpp:196
Mesquite::QualityMetric::set_metric_type
void set_metric_type(MetricType t)
This function should be used in the constructor of every concrete quality metric. ...
Definition: includeLinks/QualityMetric.hpp:291
Mesquite::MSQ_SQRT_THREE_INV
static const double MSQ_SQRT_THREE_INV
Definition: Mesquite.hpp:125
Mesquite::QualityMetric::average_metrics
double average_metrics(const double metric_values[], const int &num_values, MsqError &err)
average_metrics takes an array of length num_values and averages the contents using averaging method ...
Definition: includeLinks/QualityMetric.hpp:524
Mesquite::MsqVertex
MsqVertex is the Mesquite object that stores information about the vertices in the mesh...
Definition: includeLinks/MsqVertex.hpp:53
Mesquite::QualityMetric::set_name
void set_name(msq_std::string st)
Sets the name of this metric.
Definition: includeLinks/QualityMetric.hpp:181
Mesquite::MSQ_SQRT_THREE
static const double MSQ_SQRT_THREE
Definition: Mesquite.hpp:123