Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
I_DFT Class Reference

Class containing the target corner matrices for the context based smoothing. More...

#include <I_DFT.hpp>

Inheritance diagram for I_DFT:
Collaboration diagram for I_DFT:

Public Member Functions

 I_DFT ()
 
virtual ~I_DFT ()
 virtual destructor ensures use of polymorphism during destruction More...
 
bool evaluate_element (PatchData &pd, MsqMeshEntity *e, double &m, MsqError &err)
 Evaluate the metric for an element. More...
 
bool compute_element_analytical_gradient (PatchData &pd, MsqMeshEntity *e, MsqVertex *fv[], Vector3D g[], int nfv, double &m, MsqError &err)
 Virtual function that computes the gradient of the QualityMetric analytically. The base class implementation of this function simply prints a warning and calls compute_numerical_gradient to calculate the gradient. This is used by metric which mType is ELEMENT_BASED. For parameters, see compute_element_gradient() . More...
 
bool compute_element_analytical_hessian (PatchData &pd, MsqMeshEntity *e, MsqVertex *fv[], Vector3D g[], Matrix3D h[], int nfv, double &m, MsqError &err)
 
 I_DFT ()
 
virtual ~I_DFT ()
 virtual destructor ensures use of polymorphism during destruction More...
 
bool evaluate_element (PatchData &pd, MsqMeshEntity *e, double &m, MsqError &err)
 Evaluate the metric for an element. More...
 
bool compute_element_analytical_gradient (PatchData &pd, MsqMeshEntity *e, MsqVertex *fv[], Vector3D g[], int nfv, double &m, MsqError &err)
 Virtual function that computes the gradient of the QualityMetric analytically. The base class implementation of this function simply prints a warning and calls compute_numerical_gradient to calculate the gradient. This is used by metric which mType is ELEMENT_BASED. For parameters, see compute_element_gradient() . More...
 
bool compute_element_analytical_hessian (PatchData &pd, MsqMeshEntity *e, MsqVertex *fv[], Vector3D g[], Matrix3D h[], int nfv, double &m, MsqError &err)
 
- Public Member Functions inherited from DistanceFromTarget
virtual ~DistanceFromTarget ()
 virtual destructor ensures use of polymorphism during destruction More...
 
virtual ~DistanceFromTarget ()
 virtual destructor ensures use of polymorphism during destruction More...
 
- Public Member Functions inherited from QualityMetric
virtual ~QualityMetric ()
 
MetricType get_metric_type ()
 
void set_element_evaluation_mode (ElementEvaluationMode mode, MsqError &err)
 Sets the evaluation mode for the ELEMENT_BASED metrics. More...
 
ElementEvaluationMode get_element_evaluation_mode ()
 Returns the evaluation mode for the metric. More...
 
void set_averaging_method (AveragingMethod method, MsqError &err)
 
void set_feasible_constraint (int alpha)
 
int get_feasible_constraint ()
 Returns the feasible flag for this metric. More...
 
void set_name (msq_std::string st)
 Sets the name of this metric. More...
 
msq_std::string get_name ()
 Returns the name of this metric (as a string). More...
 
double vertex_barrier_function (double det, double delta)
 Escobar Barrier Function for Shape and Other Metrics. More...
 
virtual bool evaluate_vertex (PatchData &, MsqVertex *, double &, MsqError &err)
 Evaluate the metric for a vertex. More...
 
void set_gradient_type (GRADIENT_TYPE grad)
 Sets gradType for this metric. More...
 
void set_hessian_type (HESSIAN_TYPE ht)
 Sets hessianType for this metric. More...
 
bool compute_vertex_gradient (PatchData &pd, MsqVertex &vertex, MsqVertex *vertices[], Vector3D grad_vec[], int num_vtx, double &metric_value, MsqError &err)
 Calls compute_vertex_numerical_gradient if gradType equals NUMERCIAL_GRADIENT. Calls compute_vertex_analytical_gradient if gradType equals ANALYTICAL_GRADIENT;. More...
 
bool compute_element_gradient (PatchData &pd, MsqMeshEntity *element, MsqVertex *free_vtces[], Vector3D grad_vec[], int num_free_vtx, double &metric_value, MsqError &err)
 For MetricType == ELEMENT_BASED. Calls either compute_element_numerical_gradient() or compute_element_analytical_gradient() for gradType equal NUMERICAL_GRADIENT or ANALYTICAL_GRADIENT, respectively. More...
 
bool compute_element_gradient_expanded (PatchData &pd, MsqMeshEntity *element, MsqVertex *free_vtces[], Vector3D grad_vec[], int num_free_vtx, double &metric_value, MsqError &err)
 
bool compute_element_hessian (PatchData &pd, MsqMeshEntity *element, MsqVertex *free_vtces[], Vector3D grad_vec[], Matrix3D hessian[], int num_free_vtx, double &metric_value, MsqError &err)
 For MetricType == ELEMENT_BASED. Calls either compute_element_numerical_hessian() or compute_element_analytical_hessian() for hessianType equal NUMERICAL_HESSIAN or ANALYTICAL_HESSIAN, respectively. More...
 
void set_negate_flag (int neg)
 
int get_negate_flag ()
 Returns negateFlag. More...
 
virtual void change_metric_type (MetricType t, MsqError &err)
 
virtual ~QualityMetric ()
 
MetricType get_metric_type ()
 
void set_element_evaluation_mode (ElementEvaluationMode mode, MsqError &err)
 Sets the evaluation mode for the ELEMENT_BASED metrics. More...
 
ElementEvaluationMode get_element_evaluation_mode ()
 Returns the evaluation mode for the metric. More...
 
void set_averaging_method (AveragingMethod method, MsqError &err)
 
void set_feasible_constraint (int alpha)
 
int get_feasible_constraint ()
 Returns the feasible flag for this metric. More...
 
void set_name (msq_std::string st)
 Sets the name of this metric. More...
 
msq_std::string get_name ()
 Returns the name of this metric (as a string). More...
 
double vertex_barrier_function (double det, double delta)
 Escobar Barrier Function for Shape and Other Metrics. More...
 
virtual bool evaluate_vertex (PatchData &, MsqVertex *, double &, MsqError &err)
 Evaluate the metric for a vertex. More...
 
void set_gradient_type (GRADIENT_TYPE grad)
 Sets gradType for this metric. More...
 
void set_hessian_type (HESSIAN_TYPE ht)
 Sets hessianType for this metric. More...
 
bool compute_vertex_gradient (PatchData &pd, MsqVertex &vertex, MsqVertex *vertices[], Vector3D grad_vec[], int num_vtx, double &metric_value, MsqError &err)
 
bool compute_element_gradient (PatchData &pd, MsqMeshEntity *element, MsqVertex *free_vtces[], Vector3D grad_vec[], int num_free_vtx, double &metric_value, MsqError &err)
 For MetricType == ELEMENT_BASED. Calls either compute_element_numerical_gradient() or compute_element_analytical_gradient() for gradType equal NUMERICAL_GRADIENT or ANALYTICAL_GRADIENT, respectively. More...
 
bool compute_element_gradient_expanded (PatchData &pd, MsqMeshEntity *element, MsqVertex *free_vtces[], Vector3D grad_vec[], int num_free_vtx, double &metric_value, MsqError &err)
 
bool compute_element_hessian (PatchData &pd, MsqMeshEntity *element, MsqVertex *free_vtces[], Vector3D grad_vec[], Matrix3D hessian[], int num_free_vtx, double &metric_value, MsqError &err)
 For MetricType == ELEMENT_BASED. Calls either compute_element_numerical_hessian() or compute_element_analytical_hessian() for hessianType equal NUMERICAL_HESSIAN or ANALYTICAL_HESSIAN, respectively. More...
 
void set_negate_flag (int neg)
 
int get_negate_flag ()
 Returns negateFlag. More...
 
virtual void change_metric_type (MetricType t, MsqError &err)
 

Protected Member Functions

void p_set_alpha (double alpha)
 access function to set mAlpha More...
 
double p_get_alpha ()
 access function to get mAlpha More...
 
void p_set_beta (double beta)
 access function to set mBeta More...
 
double p_get_beta ()
 access function to get mBeta More...
 
void p_set_gamma (double gamma)
 access function to set mGamma More...
 
double p_get_gamma ()
 access function to get mGamma More...
 
void p_set_use_barrier_delta (bool use_delta)
 access function to set useBarrierDelta More...
 
bool p_get_use_barrier_delta ()
 access function to get useBarrierDelta More...
 
void p_set_alpha (double alpha)
 access function to set mAlpha More...
 
double p_get_alpha ()
 access function to get mAlpha More...
 
void p_set_beta (double beta)
 access function to set mBeta More...
 
double p_get_beta ()
 access function to get mBeta More...
 
void p_set_gamma (double gamma)
 access function to set mGamma More...
 
double p_get_gamma ()
 access function to get mGamma More...
 
void p_set_use_barrier_delta (bool use_delta)
 access function to set useBarrierDelta More...
 
bool p_get_use_barrier_delta ()
 access function to get useBarrierDelta More...
 
- Protected Member Functions inherited from DistanceFromTarget
void compute_T_matrices (MsqMeshEntity &elem, PatchData &pd, Matrix3D T[], size_t num_T, double c_k[], MsqError &err)
 For a given element, compute each corner matrix A, and given a target corner matrix W, returns $ T=AW^{-1} $ for each corner. More...
 
bool get_barrier_function (PatchData &pd, const double &tau, double &h, MsqError &err)
 
void compute_T_matrices (MsqMeshEntity &elem, PatchData &pd, Matrix3D T[], size_t num_T, double c_k[], MsqError &err)
 For a given element, compute each corner matrix A, and given a target corner matrix W, returns $ T=AW^{-1} $ for each corner. More...
 
bool get_barrier_function (PatchData &pd, const double &tau, double &h, MsqError &err)
 
- Protected Member Functions inherited from QualityMetric
 QualityMetric ()
 
void set_metric_type (MetricType t)
 This function should be used in the constructor of every concrete quality metric. More...
 
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 data member avgMethod . More...
 
double average_metric_and_weights (double metric_values[], int num_metric_values, MsqError &err)
 Given a list of metric values, calculate the average metric valude according to the current avgMethod and write into the passed metric_values array the the value weight/count to use when averaging gradient vectors for the metric. More...
 
double weighted_average_metrics (const double coef[], const double metric_values[], const int &num_values, MsqError &err)
 takes an array of coefficients and an array of metrics (both of length num_value) and averages the contents using averaging method 'method'. More...
 
bool compute_vertex_numerical_gradient (PatchData &pd, MsqVertex &vertex, MsqVertex *vertices[], Vector3D grad_vec[], int num_vtx, double &metric_value, MsqError &err)
 
bool compute_element_numerical_gradient (PatchData &pd, MsqMeshEntity *element, MsqVertex *free_vtces[], Vector3D grad_vec[], int num_free_vtx, double &metric_value, MsqError &err)
 Non-virtual function which numerically computes the gradient of a QualityMetric of a given element for a given set of free vertices on that element. This is used by metric which mType is ELEMENT_BASED. For parameters, see compute_element_gradient() . More...
 
virtual bool compute_vertex_analytical_gradient (PatchData &pd, MsqVertex &vertex, MsqVertex *vertices[], Vector3D grad_vec[], int num_vtx, double &metric_value, MsqError &err)
 Virtual function that computes the gradient of the QualityMetric analytically. The base class implementation of this function simply prints a warning and calls compute_numerical_gradient to calculate the gradient. This is used by metric which mType is VERTEX_BASED. More...
 
bool compute_element_numerical_hessian (PatchData &pd, MsqMeshEntity *element, MsqVertex *free_vtces[], Vector3D grad_vec[], Matrix3D hessian[], int num_free_vtx, double &metric_value, MsqError &err)
 
 QualityMetric ()
 
void set_metric_type (MetricType t)
 This function should be used in the constructor of every concrete quality metric. More...
 
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 data member avgMethod . More...
 
double average_metric_and_weights (double metric_values[], int num_metric_values, MsqError &err)
 Given a list of metric values, calculate the average metric valude according to the current avgMethod and write into the passed metric_values array the the value weight/count to use when averaging gradient vectors for the metric. More...
 
double weighted_average_metrics (const double coef[], const double metric_values[], const int &num_values, MsqError &err)
 takes an array of coefficients and an array of metrics (both of length num_value) and averages the contents using averaging method 'method'. More...
 
bool compute_vertex_numerical_gradient (PatchData &pd, MsqVertex &vertex, MsqVertex *vertices[], Vector3D grad_vec[], int num_vtx, double &metric_value, MsqError &err)
 
bool compute_element_numerical_gradient (PatchData &pd, MsqMeshEntity *element, MsqVertex *free_vtces[], Vector3D grad_vec[], int num_free_vtx, double &metric_value, MsqError &err)
 Non-virtual function which numerically computes the gradient of a QualityMetric of a given element for a given set of free vertices on that element. This is used by metric which mType is ELEMENT_BASED. For parameters, see compute_element_gradient() . More...
 
virtual bool compute_vertex_analytical_gradient (PatchData &pd, MsqVertex &vertex, MsqVertex *vertices[], Vector3D grad_vec[], int num_vtx, double &metric_value, MsqError &err)
 Virtual function that computes the gradient of the QualityMetric analytically. The base class implementation of this function simply prints a warning and calls compute_numerical_gradient to calculate the gradient. This is used by metric which mType is VERTEX_BASED. More...
 
bool compute_element_numerical_hessian (PatchData &pd, MsqMeshEntity *element, MsqVertex *free_vtces[], Vector3D grad_vec[], Matrix3D hessian[], int num_free_vtx, double &metric_value, MsqError &err)
 

Private Attributes

double mAlpha
 
double mBeta
 
Exponent mGamma
 
bool useBarrierDelta
 
Vector3D mNormals [4]
 
Vector3D mCoords [4]
 
Vector3D mGrads [4]
 
Vector3D mAccGrads [8]
 
Matrix3D mHessians [10]
 
Matrix3D mR
 
Matrix3D invR
 
Matrix3D mQ
 

Additional Inherited Members

- Public Types inherited from QualityMetric
enum  MetricType {
  MT_UNDEFINED, VERTEX_BASED, ELEMENT_BASED, VERTEX_BASED_FREE_ONLY,
  MT_UNDEFINED, VERTEX_BASED, ELEMENT_BASED, VERTEX_BASED_FREE_ONLY
}
 
enum  ElementEvaluationMode {
  EEM_UNDEFINED, ELEMENT_VERTICES, LINEAR_GAUSS_POINTS, QUADRATIC_GAUSS_POINTS,
  CUBIC_GAUSS_POINTS, EEM_UNDEFINED, ELEMENT_VERTICES, LINEAR_GAUSS_POINTS,
  QUADRATIC_GAUSS_POINTS, CUBIC_GAUSS_POINTS
}
 
enum  AveragingMethod {
  NONE, LINEAR, RMS, HMS,
  MINIMUM, MAXIMUM, HARMONIC, GEOMETRIC,
  SUM, SUM_SQUARED, GENERALIZED_MEAN, STANDARD_DEVIATION,
  MAX_OVER_MIN, MAX_MINUS_MIN, SUM_OF_RATIOS_SQUARED, NONE,
  LINEAR, RMS, HMS, MINIMUM,
  MAXIMUM, HARMONIC, GEOMETRIC, SUM,
  SUM_SQUARED, GENERALIZED_MEAN, STANDARD_DEVIATION, MAX_OVER_MIN,
  MAX_MINUS_MIN, SUM_OF_RATIOS_SQUARED
}
 
enum  GRADIENT_TYPE { NUMERICAL_GRADIENT, ANALYTICAL_GRADIENT, NUMERICAL_GRADIENT, ANALYTICAL_GRADIENT }
 
enum  HESSIAN_TYPE { NUMERICAL_HESSIAN, ANALYTICAL_HESSIAN, NUMERICAL_HESSIAN, ANALYTICAL_HESSIAN }
 
enum  MetricType {
  MT_UNDEFINED, VERTEX_BASED, ELEMENT_BASED, VERTEX_BASED_FREE_ONLY,
  MT_UNDEFINED, VERTEX_BASED, ELEMENT_BASED, VERTEX_BASED_FREE_ONLY
}
 
enum  ElementEvaluationMode {
  EEM_UNDEFINED, ELEMENT_VERTICES, LINEAR_GAUSS_POINTS, QUADRATIC_GAUSS_POINTS,
  CUBIC_GAUSS_POINTS, EEM_UNDEFINED, ELEMENT_VERTICES, LINEAR_GAUSS_POINTS,
  QUADRATIC_GAUSS_POINTS, CUBIC_GAUSS_POINTS
}
 
enum  AveragingMethod {
  NONE, LINEAR, RMS, HMS,
  MINIMUM, MAXIMUM, HARMONIC, GEOMETRIC,
  SUM, SUM_SQUARED, GENERALIZED_MEAN, STANDARD_DEVIATION,
  MAX_OVER_MIN, MAX_MINUS_MIN, SUM_OF_RATIOS_SQUARED, NONE,
  LINEAR, RMS, HMS, MINIMUM,
  MAXIMUM, HARMONIC, GEOMETRIC, SUM,
  SUM_SQUARED, GENERALIZED_MEAN, STANDARD_DEVIATION, MAX_OVER_MIN,
  MAX_MINUS_MIN, SUM_OF_RATIOS_SQUARED
}
 
enum  GRADIENT_TYPE { NUMERICAL_GRADIENT, ANALYTICAL_GRADIENT, NUMERICAL_GRADIENT, ANALYTICAL_GRADIENT }
 
enum  HESSIAN_TYPE { NUMERICAL_HESSIAN, ANALYTICAL_HESSIAN, NUMERICAL_HESSIAN, ANALYTICAL_HESSIAN }
 
- Protected Attributes inherited from QualityMetric
AveragingMethod avgMethod
 
int feasible
 
msq_std::string metricName
 

Detailed Description

Class containing the target corner matrices for the context based smoothing.

The form of this metric is as follows (taken from I_DFTFamilyFunctions.hpp, see that file for more detail):

              mAlpha * || A*inv(W) - mBeta * I ||_F^2                    \n

---------------------------------------------------------------—
0.5^(mGamma)*(det(A*inv(W)) + sqrt(det(A*inv(W))^2 + 4*delta^2))^(mGamma)
The default for data members (corresponding to the variables above):

  mAlpha = 1/2.0; \n
  mBeta = 1.0; \n
  mGamma = MSQ_TWO_THIRDS; \n

delta, above, is calculated using PatchData::get_barrier_delta(MsqError &err), if useBarrierDelta == true. Otherwise, delta is zero.

Definition at line 69 of file includeLinks/I_DFT.hpp.

Constructor & Destructor Documentation

I_DFT ( )
inline

Definition at line 73 of file includeLinks/I_DFT.hpp.

References QualityMetric::ANALYTICAL_GRADIENT, QualityMetric::ANALYTICAL_HESSIAN, QualityMetric::ELEMENT_BASED, QualityMetric::LINEAR, I_DFT::mAlpha, I_DFT::mBeta, I_DFT::mGamma, Mesquite::MSQ_TWO_THIRDS, QualityMetric::set_averaging_method(), QualityMetric::set_gradient_type(), QualityMetric::set_hessian_type(), QualityMetric::set_metric_type(), QualityMetric::set_name(), and I_DFT::useBarrierDelta.

74  {
75  MsqError err;
80  set_name("I_DFT");
81  mAlpha = 1/2.0;
82  mBeta = 1.0;
84  useBarrierDelta = true;
85  }
void set_averaging_method(AveragingMethod method, MsqError &err)
static const double MSQ_TWO_THIRDS
Definition: Mesquite.hpp:129
void set_hessian_type(HESSIAN_TYPE ht)
Sets hessianType for this metric.
void set_gradient_type(GRADIENT_TYPE grad)
Sets gradType for this metric.
void set_metric_type(MetricType t)
This function should be used in the constructor of every concrete quality metric. ...
void set_name(msq_std::string st)
Sets the name of this metric.

Here is the call graph for this function:

virtual ~I_DFT ( )
inlinevirtual

virtual destructor ensures use of polymorphism during destruction

Definition at line 88 of file includeLinks/I_DFT.hpp.

89  {};
I_DFT ( )
inline

Definition at line 73 of file src/QualityMetric/DFT/I_DFT.hpp.

References QualityMetric::ANALYTICAL_GRADIENT, QualityMetric::ANALYTICAL_HESSIAN, QualityMetric::ELEMENT_BASED, QualityMetric::LINEAR, I_DFT::mAlpha, I_DFT::mBeta, I_DFT::mGamma, Mesquite::MSQ_TWO_THIRDS, QualityMetric::set_averaging_method(), QualityMetric::set_gradient_type(), QualityMetric::set_hessian_type(), QualityMetric::set_metric_type(), QualityMetric::set_name(), and I_DFT::useBarrierDelta.

74  {
75  MsqError err;
80  set_name("I_DFT");
81  mAlpha = 1/2.0;
82  mBeta = 1.0;
84  useBarrierDelta = true;
85  }
void set_averaging_method(AveragingMethod method, MsqError &err)
static const double MSQ_TWO_THIRDS
Definition: Mesquite.hpp:129
void set_hessian_type(HESSIAN_TYPE ht)
Sets hessianType for this metric.
void set_gradient_type(GRADIENT_TYPE grad)
Sets gradType for this metric.
void set_metric_type(MetricType t)
This function should be used in the constructor of every concrete quality metric. ...
void set_name(msq_std::string st)
Sets the name of this metric.

Here is the call graph for this function:

virtual ~I_DFT ( )
inlinevirtual

virtual destructor ensures use of polymorphism during destruction

Definition at line 88 of file src/QualityMetric/DFT/I_DFT.hpp.

89  {};

Member Function Documentation

bool compute_element_analytical_gradient ( PatchData pd,
MsqMeshEntity element,
MsqVertex free_vtces[],
Vector3D  grad_vec[],
int  num_free_vtx,
double &  metric_value,
MsqError err 
)
virtual

Virtual function that computes the gradient of the QualityMetric analytically. The base class implementation of this function simply prints a warning and calls compute_numerical_gradient to calculate the gradient. This is used by metric which mType is ELEMENT_BASED. For parameters, see compute_element_gradient() .

If that function is not over-riden in the concrete class, the base

Parameters description, see QualityMetric::compute_element_gradient() .

Returns
true if the element is valid, false otherwise.

Reimplemented from QualityMetric.

Definition at line 174 of file QualityMetric/DFT/I_DFT.cpp.

References MsqMeshEntity::compute_corner_normals(), Mesquite::g_gdft_2(), Mesquite::g_gdft_2_v0(), Mesquite::g_gdft_2_v1(), Mesquite::g_gdft_2_v2(), Mesquite::g_gdft_3(), Mesquite::g_gdft_3_v0(), Mesquite::g_gdft_3_v1(), Mesquite::g_gdft_3_v2(), Mesquite::g_gdft_3_v3(), PatchData::get_barrier_delta(), TargetMatrix::get_cK(), PatchData::get_element_index(), MsqMeshEntity::get_element_type(), PatchData::get_vertex_array(), MsqMeshEntity::get_vertex_index_array(), Mesquite::HEXAHEDRON, i, Mesquite::inv(), I_DFT::invR, j, Mesquite::m_gdft_2(), Mesquite::m_gdft_3(), I_DFT::mAccGrads, I_DFT::mAlpha, I_DFT::mBeta, I_DFT::mCoords, I_DFT::mGamma, I_DFT::mGrads, I_DFT::mNormals, I_DFT::mQ, I_DFT::mR, Mesquite::MSQ_3RT_2_OVER_6RT_3, MSQ_ERRZERO, Mesquite::MSQ_ONE_THIRD, MSQ_SETERR, MsqError::NOT_IMPLEMENTED, Mesquite::QR(), Mesquite::QUADRILATERAL, PatchData::targetMatrices, Mesquite::TETRAHEDRON, Mesquite::TRIANGLE, I_DFT::useBarrierDelta, and MsqMeshEntity::vertex_count().

181 {
182  // Only works with the weighted average
183 
184  MsqVertex *vertices = pd.get_vertex_array(err); MSQ_ERRZERO(err);
185  EntityTopology topo = e->get_element_type();
186 
187  const size_t nv = e->vertex_count();
188  const size_t *v_i = e->get_vertex_index_array();
189 
190  size_t idx = pd.get_element_index(e);
191  const TargetMatrix *W = pd.targetMatrices.get_element_corner_tags(&pd, idx, err );
192  MSQ_ERRZERO(err);
193 
194  // Initialize constants for the metric
195  const double delta = useBarrierDelta ? pd.get_barrier_delta(err) :
196  (mGamma ? 0 : 1);
197  //const double delta = useBarrierDelta ? pd.get_barrier_delta(err) : 0;
198  MSQ_ERRZERO(err);
199 
200  const int tetInd[4][4] = {{0, 1, 2, 3}, {1, 2, 0, 3},
201  {2, 0, 1, 3}, {3, 2, 1, 0}};
202  const int hexInd[8][4] = {{0, 1, 3, 4}, {1, 2, 0, 5},
203  {2, 3, 1, 6}, {3, 0, 2, 7},
204  {4, 7, 5, 0}, {5, 4, 6, 1},
205  {6, 5, 7, 2}, {7, 6, 4, 3}};
206 
207  // Variables used for computing the metric
208  double mMetric; // Metric value
209  bool mValid; // Validity of the metric
210  int i, j, mVert;
211 
212  m = 0.0;
213  switch(topo) {
214  case TRIANGLE:
215  assert(3 == nv);
216 
217  e->compute_corner_normals( mNormals, pd, err ); MSQ_ERRZERO(err);
218 
219  // The following analytic calculation only works correctly if the
220  // normal is constant. If the normal is not constant, you need
221  // to get the gradient of the normal with respect to the vertex
222  // positions to obtain the correct values.
223 
224  if (1 == nfv) {
225  // One free vertex; use the specialized code for computing the gradient.
226 
227  g[0] = 0.0;
228  for (i = 0; i < 3; ++i) {
229  mVert = -1;
230  for (j = 0; j < 3; ++j) {
231  mCoords[j] = vertices[v_i[tetInd[i][j]]];
232  if (vertices + v_i[tetInd[i][j]] == fv[0]) {
233  mVert = j;
234  }
235  }
236 
237  if (mVert >= 0) {
239 
240  QR(mQ, mR, W[i]);
241  inv(invR, mR);
242 
243  switch(mVert) {
244  case 0:
245  mValid = g_gdft_2_v0(mMetric, mGrads[0], mCoords, mNormals[i],
246  invR, mQ, mAlpha, mGamma, delta, mBeta);
247  break;
248 
249  case 1:
250  mValid = g_gdft_2_v1(mMetric, mGrads[0], mCoords, mNormals[i],
251  invR, mQ, mAlpha, mGamma, delta, mBeta);
252  break;
253 
254  default:
255  mValid = g_gdft_2_v2(mMetric, mGrads[0], mCoords, mNormals[i],
256  invR, mQ, mAlpha, mGamma, delta, mBeta);
257  break;
258  }
259 
260  if (!mValid) return false;
261  m += W[i].get_cK() * mMetric;
262  g[0] += W[i].get_cK() * mGrads[0];
263  }
264  else {
265  // For triangles, the free vertex must appear in every element.
266  // Therefore, there these accumulations should not get used.
267 
269 
270  QR(mQ, mR, W[i]);
271  inv(invR, mR);
272 
273  mValid = m_gdft_2(mMetric, mCoords, mNormals[i],
274  invR, mQ, mAlpha, mGamma, delta, mBeta);
275  if (!mValid) return false;
276  m += W[i].get_cK() * mMetric;
277  }
278  }
279 
280  m *= MSQ_ONE_THIRD;
281  g[0] *= MSQ_ONE_THIRD;
282  }
283  else {
284  for (i = 0; i < 3; ++i) {
285  mAccGrads[i] = 0.0;
286  }
287 
288  for (i = 0; i < 3; ++i) {
289  for (j = 0; j < 3; ++j) {
290  mCoords[j] = vertices[v_i[tetInd[i][j]]];
291  }
292 
294 
295  QR(mQ, mR, W[i]);
296  inv(invR, mR);
297  mValid = g_gdft_2(mMetric, mGrads, mCoords, mNormals[i], invR, mQ,
298  mAlpha, mGamma, delta, mBeta);
299 
300  if (!mValid) return false;
301  m += W[i].get_cK() * mMetric;
302  for (j = 0; j < 3; ++j) {
303  mAccGrads[tetInd[i][j]] += W[i].get_cK() * mGrads[j];
304  }
305  }
306 
307  m *= MSQ_ONE_THIRD;
308  for (i = 0; i < 3; ++i) {
310  }
311 
312  // This is not very efficient, but is one way to select correct
313  // gradients. For gradients, info is returned only for free
314  // vertices, in the order of fv[].
315 
316  for (i = 0; i < 3; ++i) {
317  for (j = 0; j < nfv; ++j) {
318  if (vertices + v_i[i] == fv[j]) {
319  g[j] = mAccGrads[i];
320  }
321  }
322  }
323  }
324  break;
325 
326  case QUADRILATERAL:
327  assert(4 == nv);
328 
329  e->compute_corner_normals( mNormals, pd, err ); MSQ_ERRZERO(err);
330 
331  // The following analytic calculation only works correctly if the
332  // normal is constant. If the normal is not constant, you need
333  // to get the gradient of the normal with respect to the vertex
334  // positions to obtain the correct values.
335 
336  if (1 == nfv) {
337  // One free vertex; use the specialized code for computing the gradient.
338 
339  g[0] = 0.0;
340  for (i = 0; i < 4; ++i) {
341  mVert = -1;
342  for (j = 0; j < 3; ++j) {
343  mCoords[j] = vertices[v_i[hexInd[i][j]]];
344  if (vertices + v_i[hexInd[i][j]] == fv[0]) {
345  mVert = j;
346  }
347  }
348 
349  if (mVert >= 0) {
350  QR(mQ, mR, W[i]);
351  inv(invR, mR);
352 
353  switch(mVert) {
354  case 0:
355  mValid = g_gdft_2_v0(mMetric, mGrads[0], mCoords, mNormals[i],
356  invR, mQ, mAlpha, mGamma, delta, mBeta);
357  break;
358 
359  case 1:
360  mValid = g_gdft_2_v1(mMetric, mGrads[0], mCoords, mNormals[i],
361  invR, mQ, mAlpha, mGamma, delta, mBeta);
362  break;
363 
364  default:
365  mValid = g_gdft_2_v2(mMetric, mGrads[0], mCoords, mNormals[i],
366  invR, mQ, mAlpha, mGamma, delta, mBeta);
367  break;
368  }
369 
370  if (!mValid) return false;
371  m += W[i].get_cK() * mMetric;
372  g[0] += W[i].get_cK() * mGrads[0];
373  }
374  else {
375  // For quadrilaterals, the free vertex only appears in three
376  // elements. Therefore, there these accumulations are needed
377  // to get the true local objective function. Note: this code
378  // can be commented out for local codes to improve performance
379  // because you are unable to change the contributions from the
380  // elements where the free vertex does not appear. (If the
381  // weight matrices change, then the code needs to be modified.)
382 
383  QR(mQ, mR, W[i]);
384  inv(invR, mR);
385 
386  mValid = m_gdft_2(mMetric, mCoords, mNormals[i],
387  invR, mQ, mAlpha, mGamma, delta, mBeta);
388  if (!mValid) return false;
389  m += W[i].get_cK() * mMetric;
390  }
391  }
392 
393  m *= 0.25;
394  g[0] *= 0.25;
395  }
396  else {
397  for (i = 0; i < 4; ++i) {
398  mAccGrads[i] = 0.0;
399  }
400 
401  for (i = 0; i < 4; ++i) {
402  for (j = 0; j < 3; ++j) {
403  mCoords[j] = vertices[v_i[hexInd[i][j]]];
404  }
405 
406  QR(mQ, mR, W[i]);
407  inv(invR, mR);
408  mValid = g_gdft_2(mMetric, mGrads, mCoords, mNormals[i], invR, mQ,
409  mAlpha, mGamma, delta, mBeta);
410 
411  if (!mValid) return false;
412  m += W[i].get_cK() * mMetric;
413  for (j = 0; j < 3; ++j) {
414  mAccGrads[hexInd[i][j]] += W[i].get_cK() * mGrads[j];
415  }
416  }
417 
418  m *= 0.25;
419  for (i = 0; i < 4; ++i) {
420  mAccGrads[i] *= 0.25;
421  }
422 
423  // This is not very efficient, but is one way to select correct gradients
424  // For gradients, info is returned only for free vertices, in the order
425  // of fv[].
426 
427  for (i = 0; i < 4; ++i) {
428  for (j = 0; j < nfv; ++j) {
429  if (vertices + v_i[i] == fv[j]) {
430  g[j] = mAccGrads[i];
431  }
432  }
433  }
434  }
435  break;
436 
437  case TETRAHEDRON:
438  assert(4 == nv);
439 
440  if (1 == nfv) {
441  // One free vertex; use the specialized code for computing the gradient.
442 
443  g[0] = 0.0;
444  for (i = 0; i < 4; ++i) {
445  mVert = -1;
446  for (j = 0; j < 4; ++j) {
447  mCoords[j] = vertices[v_i[tetInd[i][j]]];
448  if (vertices + v_i[tetInd[i][j]] == fv[0]) {
449  mVert = j;
450  }
451  }
452 
453  if (mVert >= 0) {
454  QR(mQ, mR, W[i]);
455  inv(invR, mR);
456 
457  switch(mVert) {
458  case 0:
459  mValid = g_gdft_3_v0(mMetric, mGrads[0], mCoords, invR, mQ,
460  mAlpha, mGamma, delta, mBeta);
461  break;
462 
463  case 1:
464  mValid = g_gdft_3_v1(mMetric, mGrads[0], mCoords, invR, mQ,
465  mAlpha, mGamma, delta, mBeta);
466  break;
467 
468  case 2:
469  mValid = g_gdft_3_v2(mMetric, mGrads[0], mCoords, invR, mQ,
470  mAlpha, mGamma, delta, mBeta);
471  break;
472 
473  default:
474  mValid = g_gdft_3_v3(mMetric, mGrads[0], mCoords, invR, mQ,
475  mAlpha, mGamma, delta, mBeta);
476  break;
477  }
478 
479  if (!mValid) return false;
480  m += W[i].get_cK() * mMetric;
481  g[0] += W[i].get_cK() * mGrads[0];
482  }
483  else {
484  // For tetrahedrons, the free vertex must appear in every element.
485  // Therefore, there these accumulations should not get used.
486 
487  QR(mQ, mR, W[i]);
488  inv(invR, mR);
489  mValid = m_gdft_3(mMetric, mCoords, invR, mQ,
490  mAlpha, mGamma, delta, mBeta);
491  if (!mValid) return false;
492  m += W[i].get_cK() * mMetric;
493  }
494  }
495 
496  m *= 0.25;
497  g[0] *= 0.25;
498  }
499  else {
500  for (i = 0; i < 4; ++i) {
501  mAccGrads[i] = 0.0;
502  }
503 
504  for (i = 0; i < 4; ++i) {
505  for (j = 0; j < 4; ++j) {
506  mCoords[j] = vertices[v_i[tetInd[i][j]]];
507  }
508 
509  QR(mQ, mR, W[i]);
510  inv(invR, mR);
511  mValid = g_gdft_3(mMetric, mGrads, mCoords, invR, mQ,
512  mAlpha, mGamma, delta, mBeta);
513 
514  if (!mValid) return false;
515  m += W[i].get_cK() * mMetric;
516 
517  for (j = 0; j < 4; ++j) {
518  mAccGrads[tetInd[i][j]] += W[i].get_cK() * mGrads[j];
519  }
520  }
521 
522  m *= 0.25;
523  for (i = 0; i < 4; ++i) {
524  mAccGrads[i] *= 0.25;
525  }
526 
527  // This is not very efficient, but is one way to select correct
528  // gradients. For gradients, info is returned only for free vertices,
529  // in the order of fv[].
530 
531  for (i = 0; i < 4; ++i) {
532  for (j = 0; j < nfv; ++j) {
533  if (vertices + v_i[i] == fv[j]) {
534  g[j] = mAccGrads[i];
535  }
536  }
537  }
538  }
539  break;
540 
541  case HEXAHEDRON:
542  assert(8 == nv);
543 
544  if (1 == nfv) {
545  // One free vertex; use the specialized code for computing the gradient.
546 
547  g[0] = 0.0;
548  for (i = 0; i < 8; ++i) {
549  mVert = -1;
550  for (j = 0; j < 4; ++j) {
551  mCoords[j] = vertices[v_i[hexInd[i][j]]];
552  if (vertices + v_i[hexInd[i][j]] == fv[0]) {
553  mVert = j;
554  }
555  }
556 
557  if (mVert >= 0) {
558  QR(mQ, mR, W[i]);
559  inv(invR, mR);
560 
561  switch(mVert) {
562  case 0:
563  mValid = g_gdft_3_v0(mMetric, mGrads[0], mCoords, invR, mQ,
564  mAlpha, mGamma, delta, mBeta);
565  break;
566 
567  case 1:
568  mValid = g_gdft_3_v1(mMetric, mGrads[0], mCoords, invR, mQ,
569  mAlpha, mGamma, delta, mBeta);
570  break;
571 
572  case 2:
573  mValid = g_gdft_3_v2(mMetric, mGrads[0], mCoords, invR, mQ,
574  mAlpha, mGamma, delta, mBeta);
575  break;
576 
577  default:
578  mValid = g_gdft_3_v3(mMetric, mGrads[0], mCoords, invR, mQ,
579  mAlpha, mGamma, delta, mBeta);
580  break;
581  }
582 
583  if (!mValid) return false;
584  m += W[i].get_cK() * mMetric;
585  g[0] += W[i].get_cK() * mGrads[0];
586  }
587  else {
588  // For hexahedrons, the free vertex only appears in four elements.
589  // Therefore, there these accumulations are needed to get the
590  // true local objective function. Note: this code can be commented
591  // out for local codes to improve performance because you are
592  // unable to change the contributions from the elements where the
593  // free vertex does not appear. (If the weight matrices change,
594  // then the code needs to be modified.)
595 
596  QR(mQ, mR, W[i]);
597  inv(invR, mR);
598  mValid = m_gdft_3(mMetric, mCoords, invR, mQ,
599  mAlpha, mGamma, delta, mBeta);
600  if (!mValid) return false;
601  m += W[i].get_cK() * mMetric;
602  }
603  }
604 
605  m *= 0.125;
606  g[0] *= 0.125;
607  }
608  else {
609  for (i = 0; i < 8; ++i) {
610  mAccGrads[i] = 0.0;
611  }
612 
613  for (i = 0; i < 8; ++i) {
614  for (j = 0; j < 4; ++j) {
615  mCoords[j] = vertices[v_i[hexInd[i][j]]];
616  }
617 
618  QR(mQ, mR, W[i]);
619  inv(invR, mR);
620  mValid = g_gdft_3(mMetric, mGrads, mCoords, invR, mQ,
621  mAlpha, mGamma, delta, mBeta);
622 
623  if (!mValid) return false;
624  m += W[i].get_cK() * mMetric;
625 
626  for (j = 0; j < 4; ++j) {
627  mAccGrads[hexInd[i][j]] += W[i].get_cK() * mGrads[j];
628  }
629  }
630 
631  m *= 0.125;
632  for (i = 0; i < 8; ++i) {
633  mAccGrads[i] *= 0.125;
634  }
635 
636  // This is not very efficient, but is one way to select correct
637  // gradients. For gradients, info is returned only for free
638  // vertices, in the order of fv[].
639 
640  for (i = 0; i < 8; ++i) {
641  for (j = 0; j < nfv; ++j) {
642  if (vertices + v_i[i] == fv[j]) {
643  g[j] = mAccGrads[i];
644  }
645  }
646  }
647  }
648  break;
649 
650  default:
651  MSQ_SETERR(err)("element type not implemented.",MsqError::NOT_IMPLEMENTED);
652  return false;
653  }
654 
655  return true;
656 }
#define MSQ_ERRZERO(err)
Return zero/NULL on error.
static const double MSQ_ONE_THIRD
Definition: Mesquite.hpp:128
CornerTag< TargetMatrix > targetMatrices
Target matrix data.
EntityTopology
Definition: Mesquite.hpp:92
bool g_gdft_2_v1(double &obj, Vector3D &g_obj, const Vector3D x[3], const Vector3D &n, const Matrix3D &invR, const Matrix3D &Q, const double alpha=0.5, const Exponent &gamma=1.0, const double delta=0.0, const double beta=0.0)
bool m_gdft_3(double &obj, const Vector3D x[4], const Matrix3D &invR, const Matrix3D &Q, const double alpha=1.0/3.0, const Exponent &gamma=2.0/3.0, const double delta=0.0, const double beta=0.0)
requested functionality is not (yet) implemented
bool g_gdft_3_v2(double &obj, Vector3D &g_obj, const Vector3D x[4], const Matrix3D &invR, const Matrix3D &Q, const double alpha=1.0/3.0, const Exponent &gamma=2.0/3.0, const double delta=0.0, const double beta=0.0)
bool g_gdft_3_v0(double &obj, Vector3D &g_obj, const Vector3D x[4], const Matrix3D &invR, const Matrix3D &Q, const double alpha=1.0/3.0, const Exponent &gamma=2.0/3.0, const double delta=0.0, const double beta=0.0)
bool g_gdft_3_v3(double &obj, Vector3D &g_obj, const Vector3D x[4], const Matrix3D &invR, const Matrix3D &Q, const double alpha=1.0/3.0, const Exponent &gamma=2.0/3.0, const double delta=0.0, const double beta=0.0)
void QR(Matrix3D &Q, Matrix3D &R, const Matrix3D &A)
double get_barrier_delta(MsqError &err)
Returns delta based on the minimum and maximum corner determinant over all elements in the patch This...
size_t get_element_index(MsqMeshEntity *element)
bool g_gdft_2_v2(double &obj, Vector3D &g_obj, const Vector3D x[3], const Vector3D &n, const Matrix3D &invR, const Matrix3D &Q, const double alpha=0.5, const Exponent &gamma=1.0, const double delta=0.0, const double beta=0.0)
#define MSQ_SETERR(err)
Macro to set error - use err.clear() to clear.
bool m_gdft_2(double &obj, const Vector3D x[3], const Vector3D &n, const Matrix3D &invR, const Matrix3D &Q, const double alpha=0.5, const Exponent &gamma=1.0, const double delta=0.0, const double beta=0.0)
blockLoc i
Definition: read.cpp:79
const MsqVertex * get_vertex_array(MsqError &err) const
Returns a pointer to the start of the vertex array.
j indices j
Definition: Indexing.h:6
bool g_gdft_3_v1(double &obj, Vector3D &g_obj, const Vector3D x[4], const Matrix3D &invR, const Matrix3D &Q, const double alpha=1.0/3.0, const Exponent &gamma=2.0/3.0, const double delta=0.0, const double beta=0.0)
Class containing the target corner matrices for the context based smoothing.
void inv(Matrix3D &Ainv, const Matrix3D &A)
bool g_gdft_2_v0(double &obj, Vector3D &g_obj, const Vector3D x[3], const Vector3D &n, const Matrix3D &invR, const Matrix3D &Q, const double alpha=0.5, const Exponent &gamma=1.0, const double delta=0.0, const double beta=0.0)
MsqVertex is the Mesquite object that stores information about the vertices in the mesh...
static const double MSQ_3RT_2_OVER_6RT_3
Definition: Mesquite.hpp:130
bool g_gdft_2(double &obj, Vector3D g_obj[3], const Vector3D x[3], const Vector3D &n, const Matrix3D &invR, const Matrix3D &Q, const double alpha=0.5, const Exponent &gamma=1.0, const double delta=0.0, const double beta=0.0)
bool g_gdft_3(double &obj, Vector3D g_obj[4], const Vector3D x[4], const Matrix3D &invR, const Matrix3D &Q, const double alpha=1.0/3.0, const Exponent &gamma=2.0/3.0, const double delta=0.0, const double beta=0.0)

Here is the call graph for this function:

bool compute_element_analytical_gradient ( PatchData pd,
MsqMeshEntity element,
MsqVertex free_vtces[],
Vector3D  grad_vec[],
int  num_free_vtx,
double &  metric_value,
MsqError err 
)
virtual

Virtual function that computes the gradient of the QualityMetric analytically. The base class implementation of this function simply prints a warning and calls compute_numerical_gradient to calculate the gradient. This is used by metric which mType is ELEMENT_BASED. For parameters, see compute_element_gradient() .

If that function is not over-riden in the concrete class, the base

Parameters description, see QualityMetric::compute_element_gradient() .

Returns
true if the element is valid, false otherwise.

Reimplemented from QualityMetric.

bool compute_element_analytical_hessian ( PatchData pd,
MsqMeshEntity element,
MsqVertex free_vtces[],
Vector3D  grad_vec[],
Matrix3D  hessian[],
int  num_free_vtx,
double &  metric_value,
MsqError err 
)
virtual

If that function is not over-riden in the concrete class, the base class function makes it default to a numerical hessian.

For parameters description, see QualityMetric::compute_element_hessian() .

Returns
true if the element is valid, false otherwise.

Reimplemented from QualityMetric.

Definition at line 658 of file QualityMetric/DFT/I_DFT.cpp.

References MsqMeshEntity::compute_corner_normals(), PatchData::get_barrier_delta(), TargetMatrix::get_cK(), PatchData::get_element_index(), MsqMeshEntity::get_element_type(), PatchData::get_vertex_array(), MsqMeshEntity::get_vertex_index_array(), Mesquite::h_gdft_2(), Mesquite::h_gdft_2_v0(), Mesquite::h_gdft_2_v1(), Mesquite::h_gdft_2_v2(), Mesquite::h_gdft_3(), Mesquite::h_gdft_3_v0(), Mesquite::h_gdft_3_v1(), Mesquite::h_gdft_3_v2(), Mesquite::h_gdft_3_v3(), Mesquite::HEXAHEDRON, i, Mesquite::inv(), I_DFT::invR, j, k, Mesquite::m_gdft_2(), Mesquite::m_gdft_3(), I_DFT::mAlpha, I_DFT::mBeta, I_DFT::mCoords, I_DFT::mGamma, I_DFT::mGrads, I_DFT::mHessians, I_DFT::mNormals, I_DFT::mQ, I_DFT::mR, Mesquite::MSQ_3RT_2_OVER_6RT_3, MSQ_ERRZERO, Mesquite::MSQ_ONE_THIRD, MSQ_SETERR, MsqError::NOT_IMPLEMENTED, Mesquite::QR(), Mesquite::QUADRILATERAL, PatchData::targetMatrices, Mesquite::TETRAHEDRON, Mesquite::transpose(), Mesquite::TRIANGLE, I_DFT::useBarrierDelta, MsqMeshEntity::vertex_count(), and Matrix3D::zero().

666 {
667  // Only works with the weighted average
668 
669  MsqVertex *vertices = pd.get_vertex_array(err); MSQ_ERRZERO(err);
670  EntityTopology topo = e->get_element_type();
671 
672  const size_t nv = e->vertex_count();
673  const size_t *v_i = e->get_vertex_index_array();
674 
675  size_t idx = pd.get_element_index(e);
676  const TargetMatrix *W = pd.targetMatrices.get_element_corner_tags(&pd, idx, err );
677  MSQ_ERRZERO(err);
678 
679  // Initialize constants for the metric
680  const double delta = useBarrierDelta ? pd.get_barrier_delta(err) :
681  (mGamma ? 0 : 1);
682  //const double delta = useBarrierDelta ? pd.get_barrier_delta(err) : 0;
683  MSQ_ERRZERO(err);
684 
685  const int tetInd[4][4] = {{0, 1, 2, 3}, {1, 2, 0, 3},
686  {2, 0, 1, 3}, {3, 2, 1, 0}};
687  const int hexInd[8][4] = {{0, 1, 3, 4}, {1, 2, 0, 5},
688  {2, 3, 1, 6}, {3, 0, 2, 7},
689  {4, 7, 5, 0}, {5, 4, 6, 1},
690  {6, 5, 7, 2}, {7, 6, 4, 3}};
691 
692  // Variables used for computing the metric
693  double mMetric; // Metric value
694  bool mValid; // Validity of the metric
695  int i, j, k, l, mVert;
696  int row, col, loc;
697 
698  m = 0.0;
699  switch(topo) {
700  case TRIANGLE:
701  assert(3 == nv);
702  e->compute_corner_normals( mNormals, pd, err ); MSQ_ERRZERO(err);
703 
704  // The following analytic calculation only works correctly if the
705  // normal is constant. If the normal is not constant, you need
706  // to get the gradient of the normal with respect to the vertex
707  // positions to obtain the correct values.
708 
709  // Zero out the hessian and gradient vector
710  for (i = 0; i < 3; ++i) {
711  g[i] = 0.0;
712  }
713 
714  for (i = 0; i < 6; ++i) {
715  h[i].zero();
716  }
717 
718  if (1 == nfv) {
719  // One free vertex; use the specialized code for computing the
720  // gradient and Hessian.
721 
722  Vector3D mG;
723  Matrix3D mH;
724 
725  mG = 0.0;
726  mH.zero();
727 
728  for (i = 0; i < 3; ++i) {
729  mVert = -1;
730  for (j = 0; j < 3; ++j) {
731  mCoords[j] = vertices[v_i[tetInd[i][j]]];
732  if (vertices + v_i[tetInd[i][j]] == fv[0]) {
733  mVert = j;
734  }
735  }
736 
737  if (mVert >= 0) {
739 
740  QR(mQ, mR, W[i]);
741  inv(invR, mR);
742 
743  switch(mVert) {
744  case 0:
745  mValid = h_gdft_2_v0(mMetric, mGrads[0], mHessians[0],
746  mCoords, mNormals[i],
747  invR, mQ, mAlpha, mGamma, delta, mBeta);
748  break;
749 
750  case 1:
751  mValid = h_gdft_2_v1(mMetric, mGrads[0], mHessians[0],
752  mCoords, mNormals[i],
753  invR, mQ, mAlpha, mGamma, delta, mBeta);
754  break;
755 
756  default:
757  mValid = h_gdft_2_v2(mMetric, mGrads[0], mHessians[0],
758  mCoords, mNormals[i],
759  invR, mQ, mAlpha, mGamma, delta, mBeta);
760  break;
761  }
762 
763  if (!mValid) return false;
764  m += W[i].get_cK() * mMetric;
765  mG += W[i].get_cK() * mGrads[0];
766  mH += W[i].get_cK() * mHessians[0];
767  }
768  else {
769  // For triangles, the free vertex must appear in every element.
770  // Therefore, there these accumulations should not get used.
771 
773 
774  QR(mQ, mR, W[i]);
775  inv(invR, mR);
776 
777  mValid = m_gdft_2(mMetric, mCoords, mNormals[i],
778  invR, mQ, mAlpha, mGamma, delta, mBeta);
779  if (!mValid) return false;
780  m += W[i].get_cK() * mMetric;
781  }
782  }
783 
784  m *= MSQ_ONE_THIRD;
785  mG *= MSQ_ONE_THIRD;
786  mH *= MSQ_ONE_THIRD;
787 
788  for (i = 0; i < 3; ++i) {
789  if (vertices + v_i[i] == fv[0]) {
790  // free vertex, see next
791  g[i] = mG;
792  switch(i) {
793  case 0:
794  h[0] = mH;
795  break;
796 
797  case 1:
798  h[3] = mH;
799  break;
800 
801  default:
802  h[5] = mH;
803  break;
804  }
805  break;
806  }
807  }
808  }
809  else {
810  // Compute the metric and sum them together
811  for (i = 0; i < 3; ++i) {
812  for (j = 0; j < 3; ++j) {
813  mCoords[j] = vertices[v_i[tetInd[i][j]]];
814  }
815 
817 
818  QR(mQ, mR, W[i]);
819  inv(invR, mR);
820  mValid = h_gdft_2(mMetric, mGrads, mHessians, mCoords, mNormals[i],
821  invR, mQ, mAlpha, mGamma, delta, mBeta);
822 
823  if (!mValid) return false;
824  m += W[i].get_cK() * mMetric;
825 
826  for (j = 0; j < 3; ++j) {
827  g[tetInd[i][j]] += W[i].get_cK() * mGrads[j];
828  }
829 
830  l = 0;
831  for (j = 0; j < 3; ++j) {
832  for (k = j; k < 3; ++k) {
833  row = tetInd[i][j];
834  col = tetInd[i][k];
835 
836  if (row <= col) {
837  loc = 3*row - (row*(row+1)/2) + col;
838  h[loc] += W[i].get_cK() * mHessians[l];
839  }
840  else {
841  loc = 3*col - (col*(col+1)/2) + row;
842  h[loc] += W[i].get_cK() * transpose(mHessians[l]);
843  }
844  ++l;
845  }
846  }
847  }
848 
849  m *= MSQ_ONE_THIRD;
850  for (i = 0; i < 3; ++i) {
851  g[i] *= MSQ_ONE_THIRD;
852  }
853 
854  for (i = 0; i < 6; ++i) {
855  h[i] *= MSQ_ONE_THIRD;
856  }
857 
858  // zero out fixed elements of g
859  j = 0;
860  for (i = 0; i < 3; ++i) {
861  if (vertices + v_i[i] == fv[j]) {
862  // if free vertex, see next
863  ++j;
864  }
865  else {
866  // else zero gradient entry and hessian entries.
867  g[i] = 0.;
868 
869  switch(i) {
870  case 0:
871  h[0].zero(); h[1].zero(); h[2].zero();
872  break;
873 
874  case 1:
875  h[1].zero(); h[3].zero(); h[4].zero();
876  break;
877 
878  case 2:
879  h[2].zero(); h[4].zero(); h[5].zero();
880  break;
881  }
882  }
883  }
884  }
885  break;
886 
887  case QUADRILATERAL:
888  assert(4 == nv);
889  e->compute_corner_normals( mNormals, pd, err ); MSQ_ERRZERO(err);
890 
891  // The following analytic calculation only works correctly if the
892  // normal is constant. If the normal is not constant, you need
893  // to get the gradient of the normal with respect to the vertex
894  // positions to obtain the correct values.
895 
896  // Zero out the hessian and gradient vector
897  for (i = 0; i < 4; ++i) {
898  g[i] = 0.0;
899  }
900 
901  for (i = 0; i < 10; ++i) {
902  h[i].zero();
903  }
904 
905  if (1 == nfv) {
906  // One free vertex; use the specialized code for computing the
907  // gradient and Hessian.
908 
909  Vector3D mG;
910  Matrix3D mH;
911 
912  mG = 0.0;
913  mH.zero();
914 
915  for (i = 0; i < 4; ++i) {
916  mVert = -1;
917  for (j = 0; j < 3; ++j) {
918  mCoords[j] = vertices[v_i[hexInd[i][j]]];
919  if (vertices + v_i[hexInd[i][j]] == fv[0]) {
920  mVert = j;
921  }
922  }
923 
924  if (mVert >= 0) {
925 
926  QR(mQ, mR, W[i]);
927  inv(invR, mR);
928 
929  switch(mVert) {
930  case 0:
931  mValid = h_gdft_2_v0(mMetric, mGrads[0], mHessians[0],
932  mCoords, mNormals[i],
933  invR, mQ, mAlpha, mGamma, delta, mBeta);
934  break;
935 
936  case 1:
937  mValid = h_gdft_2_v1(mMetric, mGrads[0], mHessians[0],
938  mCoords, mNormals[i],
939  invR, mQ, mAlpha, mGamma, delta, mBeta);
940  break;
941 
942  default:
943  mValid = h_gdft_2_v2(mMetric, mGrads[0], mHessians[0],
944  mCoords, mNormals[i],
945  invR, mQ, mAlpha, mGamma, delta, mBeta);
946  break;
947  }
948 
949  if (!mValid) return false;
950  m += W[i].get_cK() * mMetric;
951  mG += W[i].get_cK() * mGrads[0];
952  mH += W[i].get_cK() * mHessians[0];
953  }
954  else {
955  // For quadrilaterals, the free vertex only appears in three
956  // elements. Therefore, there these accumulations are needed
957  // to get the true local objective function. Note: this code
958  // can be commented out for local codes to improve performance
959  // because you are unable to change the contributions from the
960  // elements where the free vertex does not appear. (If the
961  // weight matrices change, then the code needs to be modified.)
962 
963  QR(mQ, mR, W[i]);
964  inv(invR, mR);
965 
966  mValid = m_gdft_2(mMetric, mCoords, mNormals[i],
967  invR, mQ, mAlpha, mGamma, delta, mBeta);
968  if (!mValid) return false;
969  m += W[i].get_cK() * mMetric;
970  }
971  }
972 
973  m *= 0.25;
974  mG *= 0.25;
975  mH *= 0.25;
976 
977  for (i = 0; i < 4; ++i) {
978  if (vertices + v_i[i] == fv[0]) {
979  // free vertex, see next
980  g[i] = mG;
981  switch(i) {
982  case 0:
983  h[0] = mH;
984  break;
985 
986  case 1:
987  h[4] = mH;
988  break;
989 
990  case 2:
991  h[7] = mH;
992  break;
993 
994  default:
995  h[9] = mH;
996  break;
997  }
998  break;
999  }
1000  }
1001  }
1002  else {
1003  // Compute the metric and sum them together
1004  for (i = 0; i < 4; ++i) {
1005  for (j = 0; j < 3; ++j) {
1006  mCoords[j] = vertices[v_i[hexInd[i][j]]];
1007  }
1008 
1009  QR(mQ, mR, W[i]);
1010  inv(invR, mR);
1011  mValid = h_gdft_2(mMetric, mGrads, mHessians, mCoords, mNormals[i],
1012  invR, mQ, mAlpha, mGamma, delta, mBeta);
1013 
1014  if (!mValid) return false;
1015  m += W[i].get_cK() * mMetric;
1016 
1017  for (j = 0; j < 3; ++j) {
1018  g[hexInd[i][j]] += W[i].get_cK() * mGrads[j];
1019  }
1020 
1021  l = 0;
1022  for (j = 0; j < 3; ++j) {
1023  for (k = j; k < 3; ++k) {
1024  row = hexInd[i][j];
1025  col = hexInd[i][k];
1026 
1027  if (row <= col) {
1028  loc = 4*row - (row*(row+1)/2) + col;
1029  h[loc] += W[i].get_cK() * mHessians[l];
1030  }
1031  else {
1032  loc = 4*col - (col*(col+1)/2) + row;
1033  h[loc] += W[i].get_cK() * transpose(mHessians[l]);
1034  }
1035  ++l;
1036  }
1037  }
1038  }
1039 
1040  m *= 0.25;
1041  for (i = 0; i < 4; ++i) {
1042  g[i] *= 0.25;
1043  }
1044 
1045  for (i = 0; i < 10; ++i) {
1046  h[i] *= 0.25;
1047  }
1048 
1049  // zero out fixed elements of gradient and Hessian
1050  j = 0;
1051  for (i = 0; i < 4; ++i) {
1052  if (vertices + v_i[i] == fv[j]) {
1053  // if free vertex, see next
1054  ++j;
1055  }
1056  else {
1057  // else zero gradient entry and hessian entries.
1058  g[i] = 0.;
1059 
1060  switch(i) {
1061  case 0:
1062  h[0].zero(); h[1].zero(); h[2].zero(); h[3].zero();
1063  break;
1064 
1065  case 1:
1066  h[1].zero(); h[4].zero(); h[5].zero(); h[6].zero();
1067  break;
1068 
1069  case 2:
1070  h[2].zero(); h[5].zero(); h[7].zero(); h[8].zero();
1071  break;
1072 
1073  case 3:
1074  h[3].zero(); h[6].zero(); h[8].zero(); h[9].zero();
1075  break;
1076  }
1077  }
1078  }
1079  }
1080  break;
1081 
1082  case TETRAHEDRON:
1083  assert(4 == nv);
1084 
1085  // Zero out the hessian and gradient vector
1086  for (i = 0; i < 4; ++i) {
1087  g[i] = 0.0;
1088  }
1089 
1090  for (i = 0; i < 10; ++i) {
1091  h[i].zero();
1092  }
1093 
1094  if (1 == nfv) {
1095  // One free vertex; use the specialized code for computing the
1096  // gradient and Hessian.
1097 
1098  Vector3D mG;
1099  Matrix3D mH;
1100 
1101  mG = 0.0;
1102  mH.zero();
1103 
1104  for (i = 0; i < 4; ++i) {
1105  mVert = -1;
1106  for (j = 0; j < 4; ++j) {
1107  mCoords[j] = vertices[v_i[tetInd[i][j]]];
1108  if (vertices + v_i[tetInd[i][j]] == fv[0]) {
1109  mVert = j;
1110  }
1111  }
1112 
1113  if (mVert >= 0) {
1114  QR(mQ, mR, W[i]);
1115  inv(invR, mR);
1116 
1117  switch(mVert) {
1118  case 0:
1119  mValid = h_gdft_3_v0(mMetric, mGrads[0], mHessians[0],
1120  mCoords, invR, mQ,
1121  mAlpha, mGamma, delta, mBeta);
1122  break;
1123 
1124  case 1:
1125  mValid = h_gdft_3_v1(mMetric, mGrads[0], mHessians[0],
1126  mCoords, invR, mQ,
1127  mAlpha, mGamma, delta, mBeta);
1128  break;
1129 
1130  case 2:
1131  mValid = h_gdft_3_v2(mMetric, mGrads[0], mHessians[0],
1132  mCoords, invR, mQ,
1133  mAlpha, mGamma, delta, mBeta);
1134  break;
1135 
1136  default:
1137  mValid = h_gdft_3_v3(mMetric, mGrads[0], mHessians[0],
1138  mCoords, invR, mQ,
1139  mAlpha, mGamma, delta, mBeta);
1140  break;
1141  }
1142 
1143  if (!mValid) return false;
1144  m += W[i].get_cK() * mMetric;
1145  mG += W[i].get_cK() * mGrads[0];
1146  mH += W[i].get_cK() * mHessians[0];
1147  }
1148  else {
1149  // For tetrahedrons, the free vertex must appear in every element.
1150  // Therefore, there these accumulations should not get used.
1151 
1152  QR(mQ, mR, W[i]);
1153  inv(invR, mR);
1154  mValid = m_gdft_3(mMetric, mCoords, invR, mQ,
1155  mAlpha, mGamma, delta, mBeta);
1156  if (!mValid) return false;
1157  m += W[i].get_cK() * mMetric;
1158  }
1159  }
1160 
1161  m *= 0.25;
1162  mG *= 0.25;
1163  mH *= 0.25;
1164 
1165  for (i = 0; i < 4; ++i) {
1166  if (vertices + v_i[i] == fv[0]) {
1167  // free vertex, see next
1168  g[i] = mG;
1169  switch(i) {
1170  case 0:
1171  h[0] = mH;
1172  break;
1173 
1174  case 1:
1175  h[4] = mH;
1176  break;
1177 
1178  case 2:
1179  h[7] = mH;
1180  break;
1181 
1182  default:
1183  h[9] = mH;
1184  break;
1185  }
1186  break;
1187  }
1188  }
1189  }
1190  else {
1191  // Compute the metric and sum them together
1192  for (i = 0; i < 4; ++i) {
1193  for (j = 0; j < 4; ++j) {
1194  mCoords[j] = vertices[v_i[tetInd[i][j]]];
1195  }
1196 
1197  QR(mQ, mR, W[i]);
1198  inv(invR, mR);
1199  mValid = h_gdft_3(mMetric, mGrads, mHessians, mCoords, invR, mQ,
1200  mAlpha, mGamma, delta, mBeta);
1201 
1202  if (!mValid) return false;
1203 
1204  m += W[i].get_cK() * mMetric;
1205 
1206  for (j = 0; j < 4; ++j) {
1207  g[tetInd[i][j]] += W[i].get_cK() * mGrads[j];
1208  }
1209 
1210  l = 0;
1211  for (j = 0; j < 4; ++j) {
1212  for (k = j; k < 4; ++k) {
1213  row = tetInd[i][j];
1214  col = tetInd[i][k];
1215 
1216  if (row <= col) {
1217  loc = 4*row - (row*(row+1)/2) + col;
1218  h[loc] += W[i].get_cK() * mHessians[l];
1219  }
1220  else {
1221  loc = 4*col - (col*(col+1)/2) + row;
1222  h[loc] += W[i].get_cK() * transpose(mHessians[l]);
1223  }
1224  ++l;
1225  }
1226  }
1227  }
1228 
1229  m *= 0.25;
1230  for (i = 0; i < 4; ++i) {
1231  g[i] *= 0.25;
1232  }
1233 
1234  for (i = 0; i < 10; ++i) {
1235  h[i] *= 0.25;
1236  }
1237 
1238  // zero out fixed elements of g
1239  j = 0;
1240  for (i = 0; i < 4; ++i) {
1241  if (vertices + v_i[i] == fv[j]) {
1242  // if free vertex, see next
1243  ++j;
1244  }
1245  else {
1246  // else zero gradient entry and hessian entries.
1247  g[i] = 0.;
1248 
1249  switch(i) {
1250  case 0:
1251  h[0].zero(); h[1].zero(); h[2].zero(); h[3].zero();
1252  break;
1253 
1254  case 1:
1255  h[1].zero(); h[4].zero(); h[5].zero(); h[6].zero();
1256  break;
1257 
1258  case 2:
1259  h[2].zero(); h[5].zero(); h[7].zero(); h[8].zero();
1260  break;
1261 
1262  case 3:
1263  h[3].zero(); h[6].zero(); h[8].zero(); h[9].zero();
1264  break;
1265  }
1266  }
1267  }
1268  }
1269  break;
1270 
1271  case HEXAHEDRON:
1272  assert(8 == nv);
1273 
1274  // Zero out the hessian and gradient vector
1275  for (i = 0; i < 8; ++i) {
1276  g[i] = 0.0;
1277  }
1278 
1279  for (i = 0; i < 36; ++i) {
1280  h[i].zero();
1281  }
1282 
1283  if (1 == nfv) {
1284  // One free vertex; use the specialized code for computing the
1285  // gradient and Hessian.
1286 
1287  Vector3D mG;
1288  Matrix3D mH;
1289 
1290  mG = 0.0;
1291  mH.zero();
1292 
1293  for (i = 0; i < 8; ++i) {
1294  mVert = -1;
1295  for (j = 0; j < 4; ++j) {
1296  mCoords[j] = vertices[v_i[hexInd[i][j]]];
1297  if (vertices + v_i[hexInd[i][j]] == fv[0]) {
1298  mVert = j;
1299  }
1300  }
1301 
1302  if (mVert >= 0) {
1303  QR(mQ, mR, W[i]);
1304  inv(invR, mR);
1305 
1306  switch(mVert) {
1307  case 0:
1308  mValid = h_gdft_3_v0(mMetric, mGrads[0], mHessians[0],
1309  mCoords, invR, mQ,
1310  mAlpha, mGamma, delta, mBeta);
1311  break;
1312 
1313  case 1:
1314  mValid = h_gdft_3_v1(mMetric, mGrads[0], mHessians[0],
1315  mCoords, invR, mQ,
1316  mAlpha, mGamma, delta, mBeta);
1317  break;
1318 
1319  case 2:
1320  mValid = h_gdft_3_v2(mMetric, mGrads[0], mHessians[0],
1321  mCoords, invR, mQ,
1322  mAlpha, mGamma, delta, mBeta);
1323  break;
1324 
1325  default:
1326  mValid = h_gdft_3_v3(mMetric, mGrads[0], mHessians[0],
1327  mCoords, invR, mQ,
1328  mAlpha, mGamma, delta, mBeta);
1329  break;
1330  }
1331 
1332  if (!mValid) return false;
1333  m += W[i].get_cK() * mMetric;
1334  mG += W[i].get_cK() * mGrads[0];
1335  mH += W[i].get_cK() * mHessians[0];
1336  }
1337  else {
1338  // For hexahedrons, the free vertex only appears in four elements.
1339  // Therefore, there these accumulations are needed to get the
1340  // true local objective function. Note: this code can be commented
1341  // out for local codes to improve performance because you are
1342  // unable to change the contributions from the elements where the
1343  // free vertex does not appear. (If the weight matrices change,
1344  // then the code needs to be modified.)
1345 
1346  QR(mQ, mR, W[i]);
1347  inv(invR, mR);
1348  mValid = m_gdft_3(mMetric, mCoords, invR, mQ,
1349  mAlpha, mGamma, delta, mBeta);
1350  if (!mValid) return false;
1351  m += W[i].get_cK() * mMetric;
1352  }
1353  }
1354 
1355  m *= 0.125;
1356  mG *= 0.125;
1357  mH *= 0.125;
1358 
1359  for (i = 0; i < 8; ++i) {
1360  if (vertices + v_i[i] == fv[0]) {
1361  // free vertex, see next
1362  g[i] = mG;
1363  switch(i) {
1364  case 0:
1365  h[0] = mH;
1366  break;
1367 
1368  case 1:
1369  h[8] = mH;
1370  break;
1371 
1372  case 2:
1373  h[15] = mH;
1374  break;
1375 
1376  case 3:
1377  h[21] = mH;
1378  break;
1379 
1380  case 4:
1381  h[26] = mH;
1382  break;
1383 
1384  case 5:
1385  h[30] = mH;
1386  break;
1387 
1388  case 6:
1389  h[33] = mH;
1390  break;
1391 
1392  default:
1393  h[35] = mH;
1394  break;
1395  }
1396  break;
1397  }
1398  }
1399  }
1400  else {
1401  // Compute the metric and sum them together
1402  for (i = 0; i < 8; ++i) {
1403  for (j = 0; j < 4; ++j) {
1404  mCoords[j] = vertices[v_i[hexInd[i][j]]];
1405  }
1406 
1407  QR(mQ, mR, W[i]);
1408  inv(invR, mR);
1409  mValid = h_gdft_3(mMetric, mGrads, mHessians, mCoords, invR, mQ,
1410  mAlpha, mGamma, delta, mBeta);
1411 
1412  if (!mValid) return false;
1413 
1414  m += W[i].get_cK() * mMetric;
1415 
1416  for (j = 0; j < 4; ++j) {
1417  g[hexInd[i][j]] += W[i].get_cK() * mGrads[j];
1418  }
1419 
1420  l = 0;
1421  for (j = 0; j < 4; ++j) {
1422  for (k = j; k < 4; ++k) {
1423  row = hexInd[i][j];
1424  col = hexInd[i][k];
1425 
1426  if (row <= col) {
1427  loc = 8*row - (row*(row+1)/2) + col;
1428  h[loc] += W[i].get_cK() * mHessians[l];
1429  }
1430  else {
1431  loc = 8*col - (col*(col+1)/2) + row;
1432  h[loc] += W[i].get_cK() * transpose(mHessians[l]);
1433  }
1434  ++l;
1435  }
1436  }
1437  }
1438 
1439  m *= 0.125;
1440  for (i = 0; i < 8; ++i) {
1441  g[i] *= 0.125;
1442  }
1443 
1444  for (i = 0; i < 36; ++i) {
1445  h[i] *= 0.125;
1446  }
1447 
1448  // zero out fixed elements of gradient and Hessian
1449  j = 0;
1450  for (i = 0; i < 8; ++i) {
1451  if (vertices + v_i[i] == fv[j]) {
1452  // if free vertex, see next
1453  ++j;
1454  }
1455  else {
1456  // else zero gradient entry and hessian entries.
1457  g[i] = 0.;
1458 
1459  switch(i) {
1460  case 0:
1461  h[0].zero(); h[1].zero(); h[2].zero(); h[3].zero();
1462  h[4].zero(); h[5].zero(); h[6].zero(); h[7].zero();
1463  break;
1464 
1465  case 1:
1466  h[1].zero(); h[8].zero(); h[9].zero(); h[10].zero();
1467  h[11].zero(); h[12].zero(); h[13].zero(); h[14].zero();
1468  break;
1469 
1470  case 2:
1471  h[2].zero(); h[9].zero(); h[15].zero(); h[16].zero();
1472  h[17].zero(); h[18].zero(); h[19].zero(); h[20].zero();
1473  break;
1474 
1475  case 3:
1476  h[3].zero(); h[10].zero(); h[16].zero(); h[21].zero();
1477  h[22].zero(); h[23].zero(); h[24].zero(); h[25].zero();
1478  break;
1479 
1480  case 4:
1481  h[4].zero(); h[11].zero(); h[17].zero(); h[22].zero();
1482  h[26].zero(); h[27].zero(); h[28].zero(); h[29].zero();
1483  break;
1484 
1485  case 5:
1486  h[5].zero(); h[12].zero(); h[18].zero(); h[23].zero();
1487  h[27].zero(); h[30].zero(); h[31].zero(); h[32].zero();
1488  break;
1489 
1490  case 6:
1491  h[6].zero(); h[13].zero(); h[19].zero(); h[24].zero();
1492  h[28].zero(); h[31].zero(); h[33].zero(); h[34].zero();
1493  break;
1494 
1495  case 7:
1496  h[7].zero(); h[14].zero(); h[20].zero(); h[25].zero();
1497  h[29].zero(); h[32].zero(); h[34].zero(); h[35].zero();
1498  break;
1499  }
1500  }
1501  }
1502  }
1503  break;
1504 
1505  default:
1506  MSQ_SETERR(err)("element type not implemented.",MsqError::NOT_IMPLEMENTED);
1507  return false;
1508  }
1509  return true;
1510 }
#define MSQ_ERRZERO(err)
Return zero/NULL on error.
static const double MSQ_ONE_THIRD
Definition: Mesquite.hpp:128
void zero()
Sets all entries to zero (more efficient than assignement).
bool h_gdft_3(double &obj, Vector3D g_obj[4], Matrix3D h_obj[10], const Vector3D x[4], const Matrix3D &invR, const Matrix3D &Q, const double alpha=1.0/3.0, const Exponent &gamma=2.0/3.0, const double delta=0.0, const double beta=0.0)
bool h_gdft_3_v0(double &obj, Vector3D &g_obj, Matrix3D &h_obj, const Vector3D x[4], const Matrix3D &invR, const Matrix3D &Q, const double alpha=1.0/3.0, const Exponent &gamma=2.0/3.0, const double delta=0.0, const double beta=0.0)
CornerTag< TargetMatrix > targetMatrices
Target matrix data.
j indices k indices k
Definition: Indexing.h:6
bool h_gdft_2_v1(double &obj, Vector3D &g_obj, Matrix3D &h_obj, const Vector3D x[3], const Vector3D &n, const Matrix3D &invR, const Matrix3D &Q, const double alpha=0.5, const Exponent &gamma=1.0, const double delta=0.0, const double beta=0.0)
EntityTopology
Definition: Mesquite.hpp:92
Matrix3D transpose(const Matrix3D &A)
bool m_gdft_3(double &obj, const Vector3D x[4], const Matrix3D &invR, const Matrix3D &Q, const double alpha=1.0/3.0, const Exponent &gamma=2.0/3.0, const double delta=0.0, const double beta=0.0)
requested functionality is not (yet) implemented
Vector3D is the object that effeciently stores information about about three-deminsional vectors...
bool h_gdft_2_v0(double &obj, Vector3D &g_obj, Matrix3D &h_obj, const Vector3D x[3], const Vector3D &n, const Matrix3D &invR, const Matrix3D &Q, const double alpha=0.5, const Exponent &gamma=1.0, const double delta=0.0, const double beta=0.0)
bool h_gdft_2(double &obj, Vector3D g_obj[3], Matrix3D h_obj[6], const Vector3D x[3], const Vector3D &n, const Matrix3D &invR, const Matrix3D &Q, const double alpha=0.5, const Exponent &gamma=1.0, const double delta=0.0, const double beta=0.0)
void QR(Matrix3D &Q, Matrix3D &R, const Matrix3D &A)
double get_barrier_delta(MsqError &err)
Returns delta based on the minimum and maximum corner determinant over all elements in the patch This...
size_t get_element_index(MsqMeshEntity *element)
bool h_gdft_3_v2(double &obj, Vector3D &g_obj, Matrix3D &h_obj, const Vector3D x[4], const Matrix3D &invR, const Matrix3D &Q, const double alpha=1.0/3.0, const Exponent &gamma=2.0/3.0, const double delta=0.0, const double beta=0.0)
3*3 Matric class, row-oriented, 0-based [i][j] indexing.
bool h_gdft_3_v1(double &obj, Vector3D &g_obj, Matrix3D &h_obj, const Vector3D x[4], const Matrix3D &invR, const Matrix3D &Q, const double alpha=1.0/3.0, const Exponent &gamma=2.0/3.0, const double delta=0.0, const double beta=0.0)
#define MSQ_SETERR(err)
Macro to set error - use err.clear() to clear.
bool m_gdft_2(double &obj, const Vector3D x[3], const Vector3D &n, const Matrix3D &invR, const Matrix3D &Q, const double alpha=0.5, const Exponent &gamma=1.0, const double delta=0.0, const double beta=0.0)
blockLoc i
Definition: read.cpp:79
const MsqVertex * get_vertex_array(MsqError &err) const
Returns a pointer to the start of the vertex array.
j indices j
Definition: Indexing.h:6
Class containing the target corner matrices for the context based smoothing.
void inv(Matrix3D &Ainv, const Matrix3D &A)
bool h_gdft_2_v2(double &obj, Vector3D &g_obj, Matrix3D &h_obj, const Vector3D x[3], const Vector3D &n, const Matrix3D &invR, const Matrix3D &Q, const double alpha=0.5, const Exponent &gamma=1.0, const double delta=0.0, const double beta=0.0)
MsqVertex is the Mesquite object that stores information about the vertices in the mesh...
bool h_gdft_3_v3(double &obj, Vector3D &g_obj, Matrix3D &h_obj, const Vector3D x[4], const Matrix3D &invR, const Matrix3D &Q, const double alpha=1.0/3.0, const Exponent &gamma=2.0/3.0, const double delta=0.0, const double beta=0.0)
static const double MSQ_3RT_2_OVER_6RT_3
Definition: Mesquite.hpp:130

Here is the call graph for this function:

bool compute_element_analytical_hessian ( PatchData pd,
MsqMeshEntity element,
MsqVertex free_vtces[],
Vector3D  grad_vec[],
Matrix3D  hessian[],
int  num_free_vtx,
double &  metric_value,
MsqError err 
)
virtual

If that function is not over-riden in the concrete class, the base class function makes it default to a numerical hessian.

For parameters description, see QualityMetric::compute_element_hessian() .

Returns
true if the element is valid, false otherwise.

Reimplemented from QualityMetric.

bool evaluate_element ( PatchData ,
MsqMeshEntity ,
double &  ,
MsqError err 
)
virtual

Evaluate the metric for an element.

Reimplemented from QualityMetric.

Definition at line 43 of file QualityMetric/DFT/I_DFT.cpp.

References MsqMeshEntity::compute_corner_normals(), PatchData::get_barrier_delta(), TargetMatrix::get_cK(), PatchData::get_element_index(), MsqMeshEntity::get_element_type(), PatchData::get_vertex_array(), MsqMeshEntity::get_vertex_index_array(), Mesquite::HEXAHEDRON, i, Mesquite::inv(), I_DFT::invR, j, Mesquite::m_gdft_2(), Mesquite::m_gdft_3(), I_DFT::mAlpha, I_DFT::mBeta, I_DFT::mCoords, I_DFT::mGamma, I_DFT::mNormals, I_DFT::mQ, I_DFT::mR, Mesquite::MSQ_3RT_2_OVER_6RT_3, MSQ_ERRZERO, Mesquite::MSQ_ONE_THIRD, MSQ_SETERR, MsqError::NOT_IMPLEMENTED, Mesquite::QR(), Mesquite::QUADRILATERAL, PatchData::targetMatrices, Mesquite::TETRAHEDRON, Mesquite::TRIANGLE, I_DFT::useBarrierDelta, and MsqMeshEntity::vertex_count().

47 {
48  // Only works with the weighted average
49 
50  MsqError mErr;
51  MsqVertex *vertices = pd.get_vertex_array(err); MSQ_ERRZERO(err);
52 
53  EntityTopology topo = e->get_element_type();
54 
55  const size_t nv = e->vertex_count();
56  const size_t *v_i = e->get_vertex_index_array();
57 
58  size_t idx = pd.get_element_index(e);
59  const TargetMatrix *W = pd.targetMatrices.get_element_corner_tags(&pd, idx, err );
60  MSQ_ERRZERO(err);
61 
62  // Initialize constants for the metric
63  const double delta = useBarrierDelta ? pd.get_barrier_delta(err) :
64  (mGamma ? 0 : 1);
65  MSQ_ERRZERO(err);
66 
67  const int tetInd[4][4] = {{0, 1, 2, 3}, {1, 2, 0, 3},
68  {2, 0, 1, 3}, {3, 2, 1, 0}};
69  const int hexInd[8][4] = {{0, 1, 3, 4}, {1, 2, 0, 5},
70  {2, 3, 1, 6}, {3, 0, 2, 7},
71  {4, 7, 5, 0}, {5, 4, 6, 1},
72  {6, 5, 7, 2}, {7, 6, 4, 3}};
73 
74  // Variables used for computing the metric
75  double mMetric; // Metric value
76  bool mValid; // Validity of the metric
77  int i, j;
78 
79  m = 0.0;
80  switch(topo) {
81  case TRIANGLE:
82  assert(3 == nv);
83 
84  e->compute_corner_normals( mNormals, pd, err ); MSQ_ERRZERO(err);
85 
86  for (i = 0; i < 3; ++i) {
87  for (j = 0; j < 3; ++j) {
88  mCoords[j] = vertices[v_i[tetInd[i][j]]];
89  }
90 
92 
93  QR(mQ, mR, W[i]);
94  inv(invR, mR);
95  mValid = m_gdft_2(mMetric, mCoords, mNormals[i], invR, mQ,
96  mAlpha, mGamma, delta, mBeta);
97  if (!mValid) return false;
98  m += W[i].get_cK() * mMetric;
99 
100  }
101 
102  m *= MSQ_ONE_THIRD;
103  break;
104 
105  case QUADRILATERAL:
106  assert(4 == nv);
107 
108  e->compute_corner_normals( mNormals, pd, err ); MSQ_ERRZERO(err);
109 
110  for (i = 0; i < 4; ++i) {
111  for (j = 0; j < 3; ++j) {
112  mCoords[j] = vertices[v_i[hexInd[i][j]]];
113  }
114 
115  QR(mQ, mR, W[i]);
116  inv(invR, mR);
117  mValid = m_gdft_2(mMetric, mCoords, mNormals[i], invR, mQ,
118  mAlpha, mGamma, delta, mBeta);
119  if (!mValid) return false;
120  m += W[i].get_cK() * mMetric;
121  }
122 
123  m *= 0.25;
124  break;
125 
126  case TETRAHEDRON:
127  assert(4 == nv);
128 
129  for (i = 0; i < 4; ++i) {
130  for (j = 0; j < 4; ++j) {
131  mCoords[j] = vertices[v_i[tetInd[i][j]]];
132  }
133 
134  QR(mQ, mR, W[i]);
135  inv(invR, mR);
136  mValid = m_gdft_3(mMetric, mCoords, invR, mQ,
137  mAlpha, mGamma, delta, mBeta);
138 
139  if (!mValid) return false;
140  m += W[i].get_cK() * mMetric;
141  }
142 
143  m *= 0.25;
144  break;
145 
146  case HEXAHEDRON:
147  assert(8 == nv);
148 
149  for (i = 0; i < 8; ++i) {
150  for (j = 0; j < 4; ++j) {
151  mCoords[j] = vertices[v_i[hexInd[i][j]]];
152  }
153 
154  QR(mQ, mR, W[i]);
155  inv(invR, mR);
156  mValid = m_gdft_3(mMetric, mCoords, invR, mQ,
157  mAlpha, mGamma, delta, mBeta);
158 
159  if (!mValid) return false;
160  m += W[i].get_cK() * mMetric;
161  }
162 
163  m *= 0.125;
164  break;
165 
166  default:
167  MSQ_SETERR(err)("element type not implemented.",MsqError::NOT_IMPLEMENTED);
168  return false;
169  }
170 
171  return true;
172 }
#define MSQ_ERRZERO(err)
Return zero/NULL on error.
static const double MSQ_ONE_THIRD
Definition: Mesquite.hpp:128
Used to hold the error state and return it to the application.
EntityTopology
Definition: Mesquite.hpp:92
bool m_gdft_3(double &obj, const Vector3D x[4], const Matrix3D &invR, const Matrix3D &Q, const double alpha=1.0/3.0, const Exponent &gamma=2.0/3.0, const double delta=0.0, const double beta=0.0)
requested functionality is not (yet) implemented
void QR(Matrix3D &Q, Matrix3D &R, const Matrix3D &A)
#define MSQ_SETERR(err)
Macro to set error - use err.clear() to clear.
bool m_gdft_2(double &obj, const Vector3D x[3], const Vector3D &n, const Matrix3D &invR, const Matrix3D &Q, const double alpha=0.5, const Exponent &gamma=1.0, const double delta=0.0, const double beta=0.0)
blockLoc i
Definition: read.cpp:79
j indices j
Definition: Indexing.h:6
Class containing the target corner matrices for the context based smoothing.
void inv(Matrix3D &Ainv, const Matrix3D &A)
MsqVertex is the Mesquite object that stores information about the vertices in the mesh...
static const double MSQ_3RT_2_OVER_6RT_3
Definition: Mesquite.hpp:130

Here is the call graph for this function:

bool evaluate_element ( PatchData ,
MsqMeshEntity ,
double &  ,
MsqError err 
)
virtual

Evaluate the metric for an element.

Reimplemented from QualityMetric.

double p_get_alpha ( )
inlineprotected

access function to get mAlpha

Definition at line 119 of file src/QualityMetric/DFT/I_DFT.hpp.

References I_DFT::mAlpha.

120  {return mAlpha;}
double p_get_alpha ( )
inlineprotected

access function to get mAlpha

Definition at line 119 of file includeLinks/I_DFT.hpp.

References I_DFT::mAlpha.

120  {return mAlpha;}
double p_get_beta ( )
inlineprotected

access function to get mBeta

Definition at line 125 of file src/QualityMetric/DFT/I_DFT.hpp.

References I_DFT::mBeta.

126  {return mBeta;}
double p_get_beta ( )
inlineprotected

access function to get mBeta

Definition at line 125 of file includeLinks/I_DFT.hpp.

References I_DFT::mBeta.

126  {return mBeta;}
double p_get_gamma ( )
inlineprotected

access function to get mGamma

Definition at line 131 of file includeLinks/I_DFT.hpp.

References I_DFT::mGamma.

132  {return mGamma;}
double p_get_gamma ( )
inlineprotected

access function to get mGamma

Definition at line 131 of file src/QualityMetric/DFT/I_DFT.hpp.

References I_DFT::mGamma.

132  {return mGamma;}
bool p_get_use_barrier_delta ( )
inlineprotected

access function to get useBarrierDelta

Definition at line 137 of file includeLinks/I_DFT.hpp.

References I_DFT::useBarrierDelta.

138  {return useBarrierDelta;}
bool p_get_use_barrier_delta ( )
inlineprotected

access function to get useBarrierDelta

Definition at line 137 of file src/QualityMetric/DFT/I_DFT.hpp.

References I_DFT::useBarrierDelta.

138  {return useBarrierDelta;}
void p_set_alpha ( double  alpha)
inlineprotected
void p_set_alpha ( double  alpha)
inlineprotected

access function to set mAlpha

Definition at line 116 of file src/QualityMetric/DFT/I_DFT.hpp.

References I_DFT::mAlpha.

117  {mAlpha = alpha;}
void p_set_beta ( double  beta)
inlineprotected

access function to set mBeta

Definition at line 122 of file src/QualityMetric/DFT/I_DFT.hpp.

References I_DFT::mBeta.

123  {mBeta = beta;}
void p_set_beta ( double  beta)
inlineprotected
void p_set_gamma ( double  gamma)
inlineprotected

access function to set mGamma

Definition at line 128 of file src/QualityMetric/DFT/I_DFT.hpp.

References I_DFT::mGamma.

129  {mGamma = gamma;}
void p_set_gamma ( double  gamma)
inlineprotected
void p_set_use_barrier_delta ( bool  use_delta)
inlineprotected
void p_set_use_barrier_delta ( bool  use_delta)
inlineprotected

access function to set useBarrierDelta

Definition at line 134 of file src/QualityMetric/DFT/I_DFT.hpp.

References I_DFT::useBarrierDelta.

135  {useBarrierDelta = use_delta;}

Member Data Documentation

Vector3D mAccGrads
private

Definition at line 151 of file includeLinks/I_DFT.hpp.

Referenced by I_DFT::compute_element_analytical_gradient().

Matrix3D mHessians
private

Definition at line 152 of file includeLinks/I_DFT.hpp.

Referenced by I_DFT::compute_element_analytical_hessian().


The documentation for this class was generated from the following files: