two additional mesh quality metrics for mesquite
This commit is contained in:
Oliver Borm (boroli)
parent
97843285d3
commit
b8f8d3b8ac
9 changed files with 4507 additions and 0 deletions
|
|
@ -0,0 +1,765 @@
|
|||
;/* *****************************************************************
|
||||
MESQUITE -- The Mesh Quality Improvement Toolkit
|
||||
|
||||
Copyright 2004 Sandia Corporation and Argonne National
|
||||
Laboratory. Under the terms of Contract DE-AC04-94AL85000
|
||||
with Sandia Corporation, the U.S. Government retains certain
|
||||
rights in this software.
|
||||
|
||||
This library is free software; you can redistribute it and/or
|
||||
modify it under the terms of the GNU Lesser General Public
|
||||
License as published by the Free Software Foundation; either
|
||||
version 2.1 of the License, or (at your option) any later version.
|
||||
|
||||
This library is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
||||
Lesser General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU Lesser General Public License
|
||||
(lgpl.txt) along with this library; if not, write to the Free Software
|
||||
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
|
||||
|
||||
diachin2@llnl.gov, djmelan@sandia.gov, mbrewer@sandia.gov,
|
||||
pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov
|
||||
|
||||
***************************************************************** */
|
||||
//
|
||||
// ORIG-DATE: 16-May-02 at 10:26:21
|
||||
// LAST-MOD: 18-Jun-04 at 11:36:07 by Thomas Leurent
|
||||
//
|
||||
/*! \file MsqMeshEntity.cpp
|
||||
|
||||
\brief All elements in Mesquite are of type MsqMeshEntity. Their associated
|
||||
functionality is implemented in this file.
|
||||
|
||||
\author Thomas Leurent
|
||||
\author Michael Brewer
|
||||
\author Darryl Melander
|
||||
\date 2002-05-16
|
||||
*/
|
||||
|
||||
#include "Mesquite.hpp"
|
||||
#include "MsqMeshEntity.hpp"
|
||||
#include "MsqVertex.hpp"
|
||||
#include "PatchData.hpp"
|
||||
|
||||
#include <vector>
|
||||
#include <ostream>
|
||||
using std::vector;
|
||||
using std::ostream;
|
||||
using std::endl;
|
||||
|
||||
namespace MESQUITE_NS {
|
||||
|
||||
//! Gets the indices of the vertices of this element.
|
||||
//! The indices are only valid in the PatchData from which
|
||||
//! this element was retrieved.
|
||||
//! The order of the vertices is the canonical order for this
|
||||
//! element's type.
|
||||
void MsqMeshEntity::get_vertex_indices(std::vector<std::size_t> &vertices) const
|
||||
{
|
||||
vertices.resize( vertex_count() );
|
||||
std::copy( vertexIndices, vertexIndices + vertex_count(), vertices.begin() );
|
||||
}
|
||||
|
||||
//! Gets the indices of the vertices of this element.
|
||||
//! The indices are only valid in the PatchData from which
|
||||
//! this element was retrieved.
|
||||
//! The order of the vertices is the canonical order for this
|
||||
//! element's type.
|
||||
//! The indices are placed appended to the end of the list.
|
||||
//! The list is not cleared before appending this entity's vertices.
|
||||
void MsqMeshEntity::append_vertex_indices(std::vector<std::size_t> &vertex_list) const
|
||||
{
|
||||
vertex_list.insert(vertex_list.end(),
|
||||
vertexIndices,
|
||||
vertexIndices + vertex_count());
|
||||
}
|
||||
|
||||
void MsqMeshEntity::get_node_indices( std::vector<std::size_t>& indices ) const
|
||||
{
|
||||
indices.resize( node_count() );
|
||||
std::copy( vertexIndices, vertexIndices + node_count(), indices.begin() );
|
||||
}
|
||||
|
||||
void MsqMeshEntity::append_node_indices( std::vector<std::size_t>& indices ) const
|
||||
{
|
||||
indices.insert( indices.end(), vertexIndices, vertexIndices + node_count() );
|
||||
}
|
||||
|
||||
|
||||
/*! The centroid of an element containing n vertices with equal masses is located at
|
||||
\f[ \b{x} = \frac{ \sum_{i=1}^{n} \b{x}_i }{ n } \f]
|
||||
where \f$ \b{x}_i ,\, i=1,...,n\f$ are the vertices coordinates.
|
||||
*/
|
||||
void MsqMeshEntity::get_centroid(Vector3D ¢roid, const PatchData &pd, MsqError &err) const
|
||||
{
|
||||
centroid = 0.0;
|
||||
const MsqVertex* vtces = pd.get_vertex_array(err); MSQ_ERRRTN(err);
|
||||
size_t nve = vertex_count();
|
||||
for (size_t i=0; i<nve; ++i)
|
||||
centroid += vtces[vertexIndices[i]];
|
||||
centroid /= nve;
|
||||
}
|
||||
|
||||
|
||||
static inline double corner_volume( const Vector3D& v0,
|
||||
const Vector3D& v1,
|
||||
const Vector3D& v2,
|
||||
const Vector3D& v3 )
|
||||
{
|
||||
return (v1 - v0) * (v2 - v0) % (v3 - v0);
|
||||
}
|
||||
|
||||
/*!
|
||||
\brief Computes the area of the given element. Returned value is
|
||||
always non-negative. If the entity passed is not a two-dimensional
|
||||
element, an error is set.*/
|
||||
double MsqMeshEntity::compute_unsigned_area(PatchData &pd, MsqError &err) {
|
||||
const MsqVertex* verts=pd.get_vertex_array(err);MSQ_ERRZERO(err);
|
||||
double tem=0.0;
|
||||
switch (mType)
|
||||
{
|
||||
|
||||
case TRIANGLE:
|
||||
tem = ((verts[vertexIndices[1]]-verts[vertexIndices[0]])*
|
||||
(verts[vertexIndices[2]]-verts[vertexIndices[0]])).length();
|
||||
return 0.5*tem;
|
||||
|
||||
case QUADRILATERAL:
|
||||
tem = ((verts[vertexIndices[1]]-verts[vertexIndices[0]])*
|
||||
(verts[vertexIndices[3]]-verts[vertexIndices[0]])).length();
|
||||
tem += ((verts[vertexIndices[3]]-verts[vertexIndices[2]])*
|
||||
(verts[vertexIndices[1]]-verts[vertexIndices[2]])).length();
|
||||
return (tem/2.0);
|
||||
|
||||
case POLYGON:
|
||||
// assume convex
|
||||
for (unsigned i = 1; i < numVertexIndices-1; ++i)
|
||||
tem += ((verts[vertexIndices[i]] - verts[vertexIndices[0]]) *
|
||||
(verts[vertexIndices[i+1]] - verts[vertexIndices[0]])).length();
|
||||
return 0.5 * tem;
|
||||
|
||||
case TETRAHEDRON:
|
||||
return 1.0/6.0 * fabs( corner_volume( verts[vertexIndices[0]],
|
||||
verts[vertexIndices[1]],
|
||||
verts[vertexIndices[2]],
|
||||
verts[vertexIndices[3]] ) );
|
||||
|
||||
case PYRAMID: {
|
||||
Vector3D m = verts[vertexIndices[0]] + verts[vertexIndices[1]]
|
||||
+ verts[vertexIndices[2]] + verts[vertexIndices[3]];
|
||||
Vector3D t1 = verts[vertexIndices[0]] - verts[vertexIndices[2]];
|
||||
Vector3D t2 = verts[vertexIndices[1]] - verts[vertexIndices[3]];
|
||||
tem = ((t1 + t2) * (t1 - t2)) % (verts[vertexIndices[4]] - 0.25 * m);
|
||||
return (1.0/12.0) * fabs(tem);
|
||||
}
|
||||
|
||||
case PRISM: {
|
||||
tem = corner_volume( verts[vertexIndices[0]],
|
||||
verts[vertexIndices[1]],
|
||||
verts[vertexIndices[2]],
|
||||
verts[vertexIndices[3]] );
|
||||
|
||||
tem += corner_volume( verts[vertexIndices[1]],
|
||||
verts[vertexIndices[2]],
|
||||
verts[vertexIndices[3]],
|
||||
verts[vertexIndices[4]] );
|
||||
|
||||
tem += corner_volume( verts[vertexIndices[2]],
|
||||
verts[vertexIndices[3]],
|
||||
verts[vertexIndices[4]],
|
||||
verts[vertexIndices[5]] );
|
||||
|
||||
return 1.0/6.0 * fabs(tem);
|
||||
}
|
||||
|
||||
case HEXAHEDRON: {
|
||||
|
||||
tem = corner_volume( verts[vertexIndices[1]],
|
||||
verts[vertexIndices[2]],
|
||||
verts[vertexIndices[0]],
|
||||
verts[vertexIndices[5]] );
|
||||
|
||||
tem += corner_volume( verts[vertexIndices[3]],
|
||||
verts[vertexIndices[0]],
|
||||
verts[vertexIndices[2]],
|
||||
verts[vertexIndices[7]] );
|
||||
|
||||
tem += corner_volume( verts[vertexIndices[4]],
|
||||
verts[vertexIndices[7]],
|
||||
verts[vertexIndices[5]],
|
||||
verts[vertexIndices[0]] );
|
||||
|
||||
tem += corner_volume( verts[vertexIndices[6]],
|
||||
verts[vertexIndices[5]],
|
||||
verts[vertexIndices[7]],
|
||||
verts[vertexIndices[2]] );
|
||||
|
||||
tem += corner_volume( verts[vertexIndices[5]],
|
||||
verts[vertexIndices[2]],
|
||||
verts[vertexIndices[0]],
|
||||
verts[vertexIndices[7]] );
|
||||
|
||||
return (1.0/6.0) * fabs(tem);
|
||||
}
|
||||
|
||||
|
||||
|
||||
default:
|
||||
MSQ_SETERR(err)("Invalid type of element passed to compute unsigned area.",
|
||||
MsqError::UNSUPPORTED_ELEMENT);
|
||||
return 0;
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
/*!
|
||||
\brief Computes the area of the given element. Returned value can be
|
||||
negative. If the entity passed is not a two-dimensional element, an
|
||||
error is set.*/
|
||||
double MsqMeshEntity::compute_signed_area(PatchData &pd, MsqError &err) {
|
||||
const MsqVertex* verts=pd.get_vertex_array(err);MSQ_ERRZERO(err);
|
||||
double tem=0.0;
|
||||
double tem2=0.0;
|
||||
Vector3D surface_normal;
|
||||
Vector3D cross_vec;
|
||||
size_t element_index=pd.get_element_index(this);
|
||||
|
||||
switch (mType)
|
||||
{
|
||||
|
||||
case TRIANGLE:
|
||||
cross_vec=((verts[vertexIndices[1]]-verts[vertexIndices[0]])*
|
||||
(verts[vertexIndices[2]]-verts[vertexIndices[0]]));
|
||||
pd.get_domain_normal_at_element(element_index,surface_normal,err);
|
||||
MSQ_ERRZERO(err);
|
||||
tem = (cross_vec.length()/2.0);
|
||||
//if normals do not point in same general direction, negate area
|
||||
if(cross_vec%surface_normal<0){
|
||||
tem *= -1;
|
||||
}
|
||||
|
||||
return tem;
|
||||
|
||||
case QUADRILATERAL:
|
||||
cross_vec=((verts[vertexIndices[1]]-verts[vertexIndices[0]])*
|
||||
(verts[vertexIndices[3]]-verts[vertexIndices[0]]));
|
||||
pd.get_domain_normal_at_element(element_index,surface_normal,err);
|
||||
MSQ_ERRZERO(err);
|
||||
tem = (cross_vec.length()/2.0);
|
||||
//if normals do not point in same general direction, negate area
|
||||
if(cross_vec%surface_normal<0){
|
||||
tem *= -1;
|
||||
}
|
||||
cross_vec=((verts[vertexIndices[3]]-verts[vertexIndices[2]])*
|
||||
(verts[vertexIndices[1]]-verts[vertexIndices[2]]));
|
||||
tem2 = (cross_vec.length()/2.0);
|
||||
//if normals do not point in same general direction, negate area
|
||||
if(cross_vec%surface_normal<0){
|
||||
tem2 *= -1;
|
||||
//test to make sure surface normal existed
|
||||
//if(surface_normal.length_squared()<.5){
|
||||
//err.set_msg("compute_signed_area called without surface_normal available.");
|
||||
//}
|
||||
}
|
||||
return (tem + tem2);
|
||||
|
||||
case PRISM: {
|
||||
tem = corner_volume( verts[vertexIndices[0]],
|
||||
verts[vertexIndices[1]],
|
||||
verts[vertexIndices[2]],
|
||||
verts[vertexIndices[3]] );
|
||||
|
||||
tem += corner_volume( verts[vertexIndices[1]],
|
||||
verts[vertexIndices[2]],
|
||||
verts[vertexIndices[3]],
|
||||
verts[vertexIndices[4]] );
|
||||
|
||||
tem += corner_volume( verts[vertexIndices[2]],
|
||||
verts[vertexIndices[3]],
|
||||
verts[vertexIndices[4]],
|
||||
verts[vertexIndices[5]] );
|
||||
|
||||
return 1.0/6.0 * tem;
|
||||
}
|
||||
default:
|
||||
MSQ_SETERR(err)("Invalid type of element passed to compute unsigned area.",
|
||||
MsqError::UNSUPPORTED_ELEMENT);
|
||||
return 0;
|
||||
};
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
/*!Appends the indices (in the vertex array) of the vertices to connected
|
||||
to vertex_array[vertex_index] to the end of the vector vert_indices.
|
||||
The connected vertices are right-hand ordered as defined by the
|
||||
entity.
|
||||
|
||||
*/
|
||||
void MsqMeshEntity::get_connected_vertices(
|
||||
std::size_t vertex_index,
|
||||
std::vector<std::size_t> &vert_indices,
|
||||
MsqError &err)
|
||||
{
|
||||
//index is set to the index in the vertexIndices corresponding
|
||||
//to vertex_index
|
||||
int index;
|
||||
for (index = vertex_count() - 1; index >= 0; --index)
|
||||
if (vertexIndices[index] == vertex_index)
|
||||
break;
|
||||
if (index < 0)
|
||||
{
|
||||
MSQ_SETERR(err)("Invalid vertex index.", MsqError::INVALID_ARG);
|
||||
return;
|
||||
}
|
||||
|
||||
unsigned n;
|
||||
const unsigned* indices = TopologyInfo::adjacent_vertices( mType, index, n );
|
||||
if (!indices)
|
||||
MSQ_SETERR(err)("Element type not available", MsqError::INVALID_ARG);
|
||||
for (unsigned i = 0; i < n; ++i)
|
||||
vert_indices.push_back( vertexIndices[indices[i]] );
|
||||
}
|
||||
|
||||
/*! Gives the normal at the surface point corner_pt ... but if not available,
|
||||
gives the normalized cross product of corner_vec1 and corner_vec2.
|
||||
*/
|
||||
/*
|
||||
void MsqMeshEntity::compute_corner_normal(size_t corner,
|
||||
Vector3D &normal,
|
||||
PatchData &pd,
|
||||
MsqError &err)
|
||||
{
|
||||
if ( get_element_type()==TRIANGLE || get_element_type()==QUADRILATERAL )
|
||||
{
|
||||
// There are two cases where we cannot get a normal from the
|
||||
// geometry that are not errors:
|
||||
// 1) There is no domain set
|
||||
// 2) The vertex is at a degenerate point on the geometry (e.g.
|
||||
// tip of a cone.)
|
||||
|
||||
bool have_normal = false;
|
||||
|
||||
// Get normal from domain
|
||||
if (pd.domain_set())
|
||||
{
|
||||
size_t index = pd.get_element_index(this);
|
||||
pd.get_domain_normal_at_corner( index, corner, normal, err );
|
||||
MSQ_ERRRTN(err);
|
||||
|
||||
double length = normal.length();
|
||||
if (length > DBL_EPSILON)
|
||||
{
|
||||
have_normal = true;
|
||||
normal /= length;
|
||||
}
|
||||
}
|
||||
|
||||
// If failed to get normal from domain, calculate
|
||||
// from adjacent vertices.
|
||||
if ( !have_normal )
|
||||
{
|
||||
const int num_corner = this->vertex_count();
|
||||
const int prev_corner = (corner + num_corner - 1) % num_corner;
|
||||
const int next_corner = (corner + 1 ) % num_corner;
|
||||
const size_t this_idx = vertexIndices[corner];
|
||||
const size_t prev_idx = vertexIndices[prev_corner];
|
||||
const size_t next_idx = vertexIndices[next_corner];
|
||||
normal = (pd.vertex_by_index(next_idx) - pd.vertex_by_index(this_idx))
|
||||
* (pd.vertex_by_index(prev_idx) - pd.vertex_by_index(this_idx));
|
||||
normal.normalize();
|
||||
}
|
||||
}
|
||||
else
|
||||
MSQ_SETERR(err)("Should only be used for faces (tri, quads, ...).",
|
||||
MsqError::INVALID_ARG);
|
||||
}
|
||||
*/
|
||||
|
||||
void MsqMeshEntity::compute_corner_normals( Vector3D normals[],
|
||||
PatchData &pd,
|
||||
MsqError &err)
|
||||
{
|
||||
EntityTopology type = get_element_type();
|
||||
if (type != TRIANGLE && type != QUADRILATERAL && type != POLYGON)
|
||||
{
|
||||
MSQ_SETERR(err)("Should only be used for faces (tri, quads, ...).",
|
||||
MsqError::INVALID_ARG);
|
||||
return;
|
||||
}
|
||||
|
||||
|
||||
// There are two cases where we cannot get a normal from the
|
||||
// geometry that are not errors:
|
||||
// 1) There is no domain set
|
||||
// 2) The vertex is at a degenerate point on the geometry (e.g.
|
||||
// tip of a cone.)
|
||||
|
||||
// Get normal from domain
|
||||
if (pd.domain_set())
|
||||
{
|
||||
size_t index = pd.get_element_index(this);
|
||||
pd.get_domain_normals_at_corners( index, normals, err );
|
||||
MSQ_ERRRTN(err);
|
||||
}
|
||||
|
||||
// Check if normals are valid (none are valid if !pd.domain_set())
|
||||
const unsigned count = vertex_count();
|
||||
for (unsigned i = 0; i < count; ++i)
|
||||
{
|
||||
// If got valid normal from domain,
|
||||
// make it a unit vector and continue.
|
||||
if (pd.domain_set())
|
||||
{
|
||||
double length = normals[i].length();
|
||||
if (length > DBL_EPSILON)
|
||||
{
|
||||
normals[i] /= length;
|
||||
continue;
|
||||
}
|
||||
}
|
||||
|
||||
const size_t prev_idx = vertexIndices[(i + count - 1) % count];
|
||||
const size_t this_idx = vertexIndices[i];
|
||||
const size_t next_idx = vertexIndices[(i + 1) % count];
|
||||
|
||||
// Calculate normal using edges adjacent to corner
|
||||
normals[i] = (pd.vertex_by_index(next_idx) - pd.vertex_by_index(this_idx))
|
||||
* (pd.vertex_by_index(prev_idx) - pd.vertex_by_index(this_idx));
|
||||
normals[i].normalize();
|
||||
}
|
||||
}
|
||||
|
||||
ostream& operator<<( ostream& stream, const MsqMeshEntity& entity )
|
||||
{
|
||||
size_t num_vert = entity.vertex_count();
|
||||
stream << "MsqMeshEntity " << &entity << " with vertices ";
|
||||
for (size_t i = 0; i < num_vert; ++i)
|
||||
stream << entity.vertexIndices[i] << " ";
|
||||
stream << endl;
|
||||
return stream;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
/*! Get a array of indices that specifies for the given vertex
|
||||
the correct matrix map. This is used by the I_DFT point
|
||||
relaxation methods in the laplacian smoothers.
|
||||
|
||||
*/
|
||||
size_t MsqMeshEntity::get_local_matrix_map_about_vertex(
|
||||
PatchData &pd, MsqVertex* vert, size_t local_map_size,
|
||||
int* local_map, MsqError &err) const
|
||||
{
|
||||
//i iterates through elem's vertices
|
||||
int i=0;
|
||||
//index is set to the index in the vertexIndices corresponding
|
||||
//to vertex_index
|
||||
int index=-1;
|
||||
int return_val=0;
|
||||
const MsqVertex* vertex_array = pd.get_vertex_array(err);
|
||||
if(err)
|
||||
return return_val;
|
||||
|
||||
switch (mType)
|
||||
{
|
||||
case TRIANGLE:
|
||||
MSQ_SETERR(err)("Requested function not yet supported for Triangles.",
|
||||
MsqError::NOT_IMPLEMENTED);
|
||||
|
||||
break;
|
||||
|
||||
case QUADRILATERAL:
|
||||
MSQ_SETERR(err)("Requested function not yet supported for Quadrilaterals.",
|
||||
MsqError::NOT_IMPLEMENTED);
|
||||
|
||||
break;
|
||||
|
||||
case TETRAHEDRON:
|
||||
if(local_map_size<4){
|
||||
MSQ_SETERR(err)("Array of incorrect length sent to function.",
|
||||
MsqError::INVALID_ARG);
|
||||
return return_val;
|
||||
}
|
||||
return_val = 4;
|
||||
while(i<4)
|
||||
{
|
||||
if(&vertex_array[vertexIndices[i]]==vert)
|
||||
{
|
||||
index=i;
|
||||
break;
|
||||
}
|
||||
++i;
|
||||
}
|
||||
switch(index){
|
||||
case(0):
|
||||
local_map[0]=0;
|
||||
local_map[1]=1;
|
||||
local_map[2]=2;
|
||||
local_map[3]=3;
|
||||
break;
|
||||
case(1):
|
||||
local_map[0]=1;
|
||||
local_map[1]=0;
|
||||
local_map[2]=3;
|
||||
local_map[3]=2;
|
||||
break;
|
||||
case(2):
|
||||
local_map[0]=2;
|
||||
local_map[1]=3;
|
||||
local_map[2]=0;
|
||||
local_map[3]=1;
|
||||
break;
|
||||
case(3):
|
||||
local_map[0]=3;
|
||||
local_map[1]=2;
|
||||
local_map[2]=1;
|
||||
local_map[3]=0;
|
||||
break;
|
||||
default:
|
||||
local_map[0]=-1;
|
||||
local_map[1]=-1;
|
||||
local_map[2]=-1;
|
||||
local_map[3]=-1;
|
||||
};
|
||||
|
||||
break;
|
||||
|
||||
case HEXAHEDRON:
|
||||
MSQ_SETERR(err)("Requested function not yet supported for Hexahedrons.",
|
||||
MsqError::NOT_IMPLEMENTED);
|
||||
|
||||
break;
|
||||
default:
|
||||
MSQ_SETERR(err)("Element type not available", MsqError::UNSUPPORTED_ELEMENT);
|
||||
break;
|
||||
}
|
||||
return return_val;
|
||||
|
||||
}
|
||||
|
||||
|
||||
void MsqMeshEntity::check_element_orientation(
|
||||
PatchData &pd, int& inverted, int& total, MsqError &err)
|
||||
{
|
||||
NodeSet all = all_nodes( err ); MSQ_ERRRTN(err);
|
||||
unsigned i;
|
||||
|
||||
if (TopologyInfo::dimension(mType) == 2) {
|
||||
if (!pd.domain_set()) {
|
||||
total = 0;
|
||||
inverted = 0;
|
||||
return;
|
||||
}
|
||||
|
||||
const MappingFunction2D* mf = pd.get_mapping_function_2D( mType );
|
||||
if (!mf) {
|
||||
check_element_orientation_corners( pd, inverted, total, err );
|
||||
return;
|
||||
}
|
||||
|
||||
NodeSet sample = mf->sample_points( all );
|
||||
total = sample.num_nodes();
|
||||
inverted = 0;
|
||||
|
||||
if (sample.have_any_corner_node()) {
|
||||
for (i = 0; i < TopologyInfo::corners(mType); ++i)
|
||||
if (sample.corner_node(i))
|
||||
inverted += inverted_jacobian_2d( pd, all, Sample(0,i), err );
|
||||
}
|
||||
if (sample.have_any_mid_edge_node()) {
|
||||
for (i = 0; i < TopologyInfo::edges(mType); ++i)
|
||||
if (sample.mid_edge_node(i))
|
||||
inverted += inverted_jacobian_2d( pd, all, Sample(0,i), err );
|
||||
}
|
||||
if (sample.have_any_mid_face_node())
|
||||
inverted += inverted_jacobian_2d( pd, all, Sample(2,0), err );
|
||||
}
|
||||
else {
|
||||
const MappingFunction3D* mf = pd.get_mapping_function_3D( mType );
|
||||
if (!mf) {
|
||||
check_element_orientation_corners( pd, inverted, total, err );
|
||||
return;
|
||||
}
|
||||
|
||||
NodeSet sample = mf->sample_points( all );
|
||||
total = sample.num_nodes();
|
||||
inverted = 0;
|
||||
|
||||
if (sample.have_any_corner_node()) {
|
||||
for (i = 0; i < TopologyInfo::corners(mType); ++i)
|
||||
if (sample.corner_node(i))
|
||||
inverted += inverted_jacobian_3d( pd, all, Sample(0,i), err );
|
||||
}
|
||||
if (sample.have_any_mid_edge_node()) {
|
||||
for (i = 0; i < TopologyInfo::edges(mType); ++i)
|
||||
if (sample.mid_edge_node(i))
|
||||
inverted += inverted_jacobian_3d( pd, all, Sample(1,i), err );
|
||||
}
|
||||
if (sample.have_any_mid_face_node()) {
|
||||
for (i = 0; i < TopologyInfo::edges(mType); ++i)
|
||||
if (sample.mid_face_node(i))
|
||||
inverted += inverted_jacobian_3d( pd, all, Sample(1,i), err );
|
||||
}
|
||||
if (sample.have_any_mid_region_node()) {
|
||||
inverted += inverted_jacobian_3d( pd, all, Sample(3,0), err );
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
bool
|
||||
MsqMeshEntity::inverted_jacobian_3d( PatchData& pd, NodeSet nodes, Sample sample, MsqError& err )
|
||||
{
|
||||
MsqMatrix<3,3> J;
|
||||
MsqVector<3> junk[27];
|
||||
size_t junk2[27], junk3;
|
||||
assert(node_count() <= 27);
|
||||
|
||||
const MappingFunction3D* mf = pd.get_mapping_function_3D( mType );
|
||||
mf->jacobian( pd, pd.get_element_index(this), nodes,
|
||||
sample, junk2, junk, junk3, J, err );
|
||||
MSQ_ERRZERO(err);
|
||||
return det(J) <= 0.0;
|
||||
}
|
||||
|
||||
bool
|
||||
MsqMeshEntity::inverted_jacobian_2d( PatchData& pd, NodeSet nodes, Sample sample, MsqError& err )
|
||||
{
|
||||
MsqMatrix<3,2> J;
|
||||
MsqVector<2> junk[9];
|
||||
size_t junk2[9], junk3;
|
||||
assert(node_count() <= 9);
|
||||
|
||||
const MappingFunction2D* mf = pd.get_mapping_function_2D( mType );
|
||||
mf->jacobian( pd, pd.get_element_index(this), nodes,
|
||||
sample, junk2, junk, junk3, J, err );
|
||||
MSQ_ERRZERO(err);
|
||||
|
||||
Vector3D norm;
|
||||
pd.get_domain_normal_at_sample( pd.get_element_index(this), sample, norm, err );
|
||||
MSQ_ERRZERO(err);
|
||||
|
||||
MsqVector<3> N(&norm[0]);
|
||||
MsqVector<3> cross = J.column(0) * J.column(1);
|
||||
return (cross % N < 0.0);
|
||||
}
|
||||
|
||||
NodeSet
|
||||
MsqMeshEntity::all_nodes( MsqError& err ) const
|
||||
{
|
||||
bool mid_edge, mid_face, mid_vol;
|
||||
TopologyInfo::higher_order( mType, node_count(), mid_edge, mid_face, mid_vol, err );
|
||||
NodeSet result;
|
||||
result.set_all_corner_nodes( mType );
|
||||
if (mid_edge)
|
||||
result.set_all_mid_edge_nodes( mType );
|
||||
if (mid_face)
|
||||
result.set_all_mid_face_nodes( mType );
|
||||
if (mid_vol)
|
||||
result.set_mid_region_node( mType );
|
||||
return result;
|
||||
}
|
||||
|
||||
void MsqMeshEntity::check_element_orientation_corners(
|
||||
PatchData &pd,
|
||||
int& inverted,
|
||||
int& total,
|
||||
MsqError &err)
|
||||
{
|
||||
total = inverted = 0;
|
||||
|
||||
if (node_count() > vertex_count()) {
|
||||
MSQ_SETERR(err)("Cannot perform inversion test for higher-order element"
|
||||
" without mapping function.", MsqError::UNSUPPORTED_ELEMENT );
|
||||
return;
|
||||
}
|
||||
|
||||
const MsqVertex *vertices = pd.get_vertex_array(err); MSQ_ERRRTN(err);
|
||||
|
||||
// Hex element descriptions
|
||||
static const int locs_hex[8][4] = {{0, 1, 3, 4},
|
||||
{1, 2, 0, 5},
|
||||
{2, 3, 1, 6},
|
||||
{3, 0, 2, 7},
|
||||
{4, 7, 5, 0},
|
||||
{5, 4, 6, 1},
|
||||
{6, 5, 7, 2},
|
||||
{7, 6, 4, 3}};
|
||||
|
||||
const Vector3D d_con(1.0, 1.0, 1.0);
|
||||
|
||||
int i;
|
||||
Vector3D coord_vectors[3];
|
||||
Vector3D center_vector;
|
||||
|
||||
switch(mType) {
|
||||
case TRIANGLE:
|
||||
|
||||
if (!pd.domain_set())
|
||||
return;
|
||||
|
||||
pd.get_domain_normal_at_element(this, coord_vectors[2], err); MSQ_ERRRTN(err);
|
||||
coord_vectors[2] = coord_vectors[2] / coord_vectors[2].length();// Need unit normal
|
||||
center_vector = vertices[vertexIndices[0]];
|
||||
coord_vectors[0] = vertices[vertexIndices[1]]-center_vector;
|
||||
coord_vectors[1] = vertices[vertexIndices[2]]-center_vector;
|
||||
total = 1;
|
||||
inverted = (coord_vectors[0]%(coord_vectors[1]*coord_vectors[2] ) <= 0.0);
|
||||
break;
|
||||
|
||||
case QUADRILATERAL:
|
||||
|
||||
if (!pd.domain_set())
|
||||
return;
|
||||
|
||||
for (i = 0; i < 4; ++i) {
|
||||
|
||||
pd.get_domain_normal_at_element(this, coord_vectors[2], err); MSQ_ERRRTN(err);
|
||||
coord_vectors[2] = coord_vectors[2] / coord_vectors[2].length();// Need unit normal
|
||||
center_vector = vertices[vertexIndices[locs_hex[i][0]]];
|
||||
coord_vectors[0] = vertices[vertexIndices[locs_hex[i][1]]]-center_vector;
|
||||
coord_vectors[1] = vertices[vertexIndices[locs_hex[i][2]]]-center_vector;
|
||||
++total;
|
||||
inverted += (coord_vectors[0]%(coord_vectors[1]*coord_vectors[2] ) <= 0.0);
|
||||
}
|
||||
break;
|
||||
|
||||
case TETRAHEDRON:
|
||||
center_vector = vertices[vertexIndices[0]];
|
||||
coord_vectors[0] = vertices[vertexIndices[1]]-center_vector;
|
||||
coord_vectors[1] = vertices[vertexIndices[2]]-center_vector;
|
||||
coord_vectors[2] = vertices[vertexIndices[3]]-center_vector;
|
||||
total = 1;
|
||||
inverted = ( coord_vectors[0]%(coord_vectors[1]*coord_vectors[2] ) <= 0.0);
|
||||
break;
|
||||
|
||||
default: // generic code for 3D elements
|
||||
{
|
||||
size_t num_corners = corner_count();
|
||||
unsigned num_adj;
|
||||
const unsigned* adj_idx;
|
||||
for (unsigned j = 0; j < num_corners; ++j)
|
||||
{
|
||||
adj_idx = TopologyInfo::adjacent_vertices( mType, j, num_adj );
|
||||
if (3 != num_adj)
|
||||
{
|
||||
MSQ_SETERR(err)("Unsupported element type.", MsqError::INVALID_ARG);
|
||||
return;
|
||||
}
|
||||
|
||||
center_vector = vertices[vertexIndices[j]];
|
||||
coord_vectors[0] = vertices[vertexIndices[adj_idx[0]]] - center_vector;
|
||||
coord_vectors[1] = vertices[vertexIndices[adj_idx[1]]] - center_vector;
|
||||
coord_vectors[2] = vertices[vertexIndices[adj_idx[2]]] - center_vector;
|
||||
++total;
|
||||
inverted += (coord_vectors[0]%(coord_vectors[1]*coord_vectors[2] ) <= 0.0);
|
||||
}
|
||||
break;
|
||||
}
|
||||
} // end switch over element type
|
||||
}
|
||||
|
||||
} // namespace Mesquite
|
||||
File diff suppressed because it is too large
Load diff
File diff suppressed because it is too large
Load diff
|
|
@ -0,0 +1,211 @@
|
|||
/* *****************************************************************
|
||||
MESQUITE -- The Mesh Quality Improvement Toolkit
|
||||
|
||||
Copyright 2004 Sandia Corporation and Argonne National
|
||||
Laboratory. Under the terms of Contract DE-AC04-94AL85000
|
||||
with Sandia Corporation, the U.S. Government retains certain
|
||||
rights in this software.
|
||||
|
||||
This library is free software; you can redistribute it and/or
|
||||
modify it under the terms of the GNU Lesser General Public
|
||||
License as published by the Free Software Foundation; either
|
||||
version 2.1 of the License, or (at your option) any later version.
|
||||
|
||||
This library is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
||||
Lesser General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU Lesser General Public License
|
||||
(lgpl.txt) along with this library; if not, write to the Free Software
|
||||
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
|
||||
|
||||
diachin2@llnl.gov, djmelan@sandia.gov, mbrewer@sandia.gov,
|
||||
pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov
|
||||
|
||||
***************************************************************** */
|
||||
/*!
|
||||
\file AspectRatioGammaQualityMetric.cpp
|
||||
|
||||
\brief Evaluates the Aspect Ratio Gamma metric for two- and
|
||||
three-diminsional simplicial elements.
|
||||
\author Michael Brewer
|
||||
\date 2002-05-09
|
||||
*/
|
||||
|
||||
#include "AspectRatioGammaQualityMetric.hpp"
|
||||
#include <math.h>
|
||||
#include "Vector3D.hpp"
|
||||
#include "MsqMeshEntity.hpp"
|
||||
#include "PatchData.hpp"
|
||||
|
||||
#include <vector>
|
||||
using std::vector;
|
||||
|
||||
using namespace Mesquite;
|
||||
|
||||
const double fourDivRootThree = 4.0/sqrt(3.0);
|
||||
const double twelveDivRootTwo = 12.0/sqrt(2.0);
|
||||
const double sixDivRootTwo = 6.0/sqrt(2.0);
|
||||
|
||||
|
||||
std::string AspectRatioGammaQualityMetric::get_name() const
|
||||
{ return "AspectRatioGamma"; }
|
||||
|
||||
int AspectRatioGammaQualityMetric::get_negate_flag() const
|
||||
{ return 1; }
|
||||
|
||||
//note that we can define this metric for other element types?
|
||||
//!Evaluate aspect ratio gamma on ``element''
|
||||
bool AspectRatioGammaQualityMetric::evaluate(PatchData &pd,
|
||||
size_t elem_index,
|
||||
double &fval,
|
||||
MsqError &err)
|
||||
{
|
||||
MsqMeshEntity& element = pd.element_by_index( elem_index );
|
||||
EntityTopology entity = element.get_element_type();
|
||||
double vol;
|
||||
Vector3D cross, normal(0,0,0);
|
||||
fval = MSQ_MAX_CAP;
|
||||
|
||||
//get element's nodes
|
||||
vert.clear();
|
||||
pd.get_element_vertex_coordinates(elem_index, vert, err); MSQ_ERRZERO(err);
|
||||
|
||||
switch(entity)
|
||||
{
|
||||
case TRIANGLE:
|
||||
//area
|
||||
cross = (vert[1] - vert[0]) * (vert[2] - vert[0]);
|
||||
vol= cross.length() / 2.0;
|
||||
if (vol < MSQ_MIN)
|
||||
return false;
|
||||
|
||||
if (pd.domain_set()) { // need domain to check for inverted elements
|
||||
pd.get_domain_normal_at_corner( elem_index, 0, normal, err ); MSQ_ERRZERO(err);
|
||||
if ((cross % normal) < -MSQ_MIN)
|
||||
return false;
|
||||
}
|
||||
|
||||
// sum squares of edge lengths
|
||||
fval = ((vert[1] - vert[0]).length_squared()
|
||||
+ (vert[2] - vert[0]).length_squared()
|
||||
+ (vert[1] - vert[2]).length_squared());
|
||||
// average
|
||||
fval /= 3.0;
|
||||
|
||||
//normalize to equil. and div by area
|
||||
fval /= vol * fourDivRootThree;
|
||||
|
||||
break;
|
||||
|
||||
case TETRAHEDRON:
|
||||
vol = (vert[1] - vert[0]) % ((vert[2] - vert[0]) * (vert[3] - vert[0])) / 6.0;
|
||||
if (vol < MSQ_MIN) // zero for degenerate and negative for inverted
|
||||
return false;
|
||||
|
||||
// sum squares of edge lengths
|
||||
fval = (vert[1] - vert[0]).length_squared()
|
||||
+ (vert[2] - vert[0]).length_squared()
|
||||
+ (vert[3] - vert[0]).length_squared()
|
||||
+ (vert[2] - vert[1]).length_squared()
|
||||
+ (vert[3] - vert[1]).length_squared()
|
||||
+ (vert[3] - vert[2]).length_squared();
|
||||
// average
|
||||
fval /= 6.0;
|
||||
|
||||
fval *= sqrt(fval);
|
||||
//normalize to equil. and div by volume
|
||||
fval /= vol * twelveDivRootTwo;
|
||||
|
||||
break;
|
||||
|
||||
case PRISM:
|
||||
// ideal shape is a prism with equilateral triangle base and equal length height
|
||||
// this volume is always positive with this method
|
||||
vol = element.compute_signed_area(pd, err);
|
||||
if (vol< MSQ_MIN) // zero for degenerate and negative for inverted
|
||||
return false;
|
||||
|
||||
fval = (vert[1] - vert[0]).length_squared()
|
||||
+ (vert[2] - vert[0]).length_squared()
|
||||
+ (vert[1] - vert[2]).length_squared()
|
||||
+ (vert[4] - vert[3]).length_squared()
|
||||
+ (vert[5] - vert[3]).length_squared()
|
||||
+ (vert[4] - vert[5]).length_squared()
|
||||
+ (vert[5] - vert[2]).length_squared()
|
||||
+ (vert[3] - vert[0]).length_squared()
|
||||
+ (vert[4] - vert[1]).length_squared();
|
||||
|
||||
// average L2 edge length
|
||||
fval /= 9.0;
|
||||
|
||||
fval *= sqrt(fval);
|
||||
//normalize to equil. and div by volume
|
||||
fval /= vol * fourDivRootThree;
|
||||
|
||||
break;
|
||||
|
||||
case HEXAHEDRON:
|
||||
// cube is an ideal shape
|
||||
// this volume is always positive with this method
|
||||
vol = element.compute_unsigned_area(pd, err);;
|
||||
if (vol< MSQ_MIN) // zero for degenerate and negative for inverted
|
||||
return false;
|
||||
|
||||
fval = (vert[1] - vert[0]).length_squared()
|
||||
+ (vert[2] - vert[1]).length_squared()
|
||||
+ (vert[3] - vert[2]).length_squared()
|
||||
+ (vert[3] - vert[0]).length_squared()
|
||||
+ (vert[5] - vert[4]).length_squared()
|
||||
+ (vert[6] - vert[5]).length_squared()
|
||||
+ (vert[7] - vert[6]).length_squared()
|
||||
+ (vert[7] - vert[4]).length_squared()
|
||||
+ (vert[6] - vert[2]).length_squared()
|
||||
+ (vert[5] - vert[1]).length_squared()
|
||||
+ (vert[4] - vert[0]).length_squared()
|
||||
+ (vert[7] - vert[3]).length_squared();
|
||||
|
||||
// average L2 edge length
|
||||
fval /= 12.0;
|
||||
|
||||
fval *= sqrt(fval);
|
||||
//normalize to equil. and div by volume
|
||||
fval /= vol;
|
||||
|
||||
break;
|
||||
|
||||
case PYRAMID:
|
||||
// Johnson Pyramid with equilateral triangle sides and square base
|
||||
// this volume is always positive with this method
|
||||
vol = element.compute_unsigned_area(pd, err);
|
||||
if (vol< MSQ_MIN) // zero for degenerate and negative for inverted
|
||||
return false;
|
||||
|
||||
fval = (vert[1] - vert[0]).length_squared()
|
||||
+ (vert[2] - vert[1]).length_squared()
|
||||
+ (vert[3] - vert[2]).length_squared()
|
||||
+ (vert[0] - vert[3]).length_squared()
|
||||
+ (vert[4] - vert[0]).length_squared()
|
||||
+ (vert[4] - vert[1]).length_squared()
|
||||
+ (vert[4] - vert[2]).length_squared()
|
||||
+ (vert[4] - vert[3]).length_squared();
|
||||
|
||||
// average L2 edge length
|
||||
fval /= 8.0;
|
||||
|
||||
fval *= sqrt(fval);
|
||||
//normalize to equil. and div by volume
|
||||
fval /= vol * sixDivRootTwo;
|
||||
|
||||
break ;
|
||||
|
||||
default:
|
||||
MSQ_SETERR(err)(MsqError::UNSUPPORTED_ELEMENT,
|
||||
"Entity type %d is not valid for Aspect Ratio Gamma\n",
|
||||
(int)entity);
|
||||
return false;
|
||||
};
|
||||
|
||||
return true;
|
||||
}
|
||||
|
|
@ -0,0 +1,21 @@
|
|||
SRC_QUALITYMETRIC_SHAPE_DIR = src/QualityMetric/Shape
|
||||
MSQ_INCLUDES += -I$(top_srcdir)/$(SRC_QUALITYMETRIC_SHAPE_DIR)
|
||||
INCLUDED_MAKEFILES += $(SRC_QUALITYMETRIC_SHAPE_DIR)/Makefile.am
|
||||
|
||||
MSQ_SRCS += \
|
||||
$(SRC_QUALITYMETRIC_SHAPE_DIR)/AspectRatioGammaQualityMetric.cpp \
|
||||
$(SRC_QUALITYMETRIC_SHAPE_DIR)/ConditionNumberQualityMetric.cpp \
|
||||
$(SRC_QUALITYMETRIC_SHAPE_DIR)/IdealWeightInverseMeanRatio.cpp \
|
||||
$(SRC_QUALITYMETRIC_SHAPE_DIR)/IdealWeightMeanRatio.cpp \
|
||||
$(SRC_QUALITYMETRIC_SHAPE_DIR)/VertexConditionNumberQualityMetric.cpp \
|
||||
$(SRC_QUALITYMETRIC_SHAPE_DIR)/VolumeRatioQualityMetric.cpp
|
||||
|
||||
MSQ_HDRS += \
|
||||
$(SRC_QUALITYMETRIC_SHAPE_DIR)/AspectRatioGammaQualityMetric.hpp \
|
||||
$(SRC_QUALITYMETRIC_SHAPE_DIR)/ConditionNumberFunctions.hpp \
|
||||
$(SRC_QUALITYMETRIC_SHAPE_DIR)/ConditionNumberQualityMetric.hpp \
|
||||
$(SRC_QUALITYMETRIC_SHAPE_DIR)/IdealWeightInverseMeanRatio.hpp \
|
||||
$(SRC_QUALITYMETRIC_SHAPE_DIR)/IdealWeightMeanRatio.hpp \
|
||||
$(SRC_QUALITYMETRIC_SHAPE_DIR)/MeanRatioFunctions.hpp \
|
||||
$(SRC_QUALITYMETRIC_SHAPE_DIR)/VertexConditionNumberQualityMetric.hpp \
|
||||
$(SRC_QUALITYMETRIC_SHAPE_DIR)/VolumeRatioQualityMetric.hpp
|
||||
|
|
@ -0,0 +1,108 @@
|
|||
/* *****************************************************************
|
||||
MESQUITE -- The Mesh Quality Improvement Toolkit
|
||||
|
||||
Copyright 2011 Oliver Borm
|
||||
|
||||
This library is free software; you can redistribute it and/or
|
||||
modify it under the terms of the GNU Lesser General Public
|
||||
License as published by the Free Software Foundation; either
|
||||
version 2.1 of the License, or (at your option) any later version.
|
||||
|
||||
This library is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
||||
Lesser General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU Lesser General Public License
|
||||
(lgpl.txt) along with this library; if not, write to the Free Software
|
||||
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
|
||||
|
||||
oli.borm@web.de
|
||||
|
||||
***************************************************************** */
|
||||
/*!
|
||||
\file VolumeRatioQualityMetric.cpp
|
||||
|
||||
\brief Evaluates the Volume Expansion Ratio metric for two- and
|
||||
three-diminsional elements.
|
||||
\author Oliver Borm
|
||||
\date 2011-01-15
|
||||
*/
|
||||
|
||||
#include "VolumeRatioQualityMetric.hpp"
|
||||
#include <math.h>
|
||||
#include "MsqMeshEntity.hpp"
|
||||
#include "PatchData.hpp"
|
||||
|
||||
#include <vector>
|
||||
using std::vector;
|
||||
|
||||
using namespace Mesquite;
|
||||
|
||||
std::string VolumeRatioQualityMetric::get_name() const
|
||||
{ return "VolumeRatio"; }
|
||||
|
||||
int VolumeRatioQualityMetric::get_negate_flag() const
|
||||
{ return -1; }
|
||||
|
||||
//!Evaluate expansion volume ratio on ``element''
|
||||
bool VolumeRatioQualityMetric::evaluate(PatchData &pd,
|
||||
size_t elem_index,
|
||||
double &fval,
|
||||
MsqError &err)
|
||||
{
|
||||
MsqMeshEntity& ownElement = pd.element_by_index( elem_index );
|
||||
|
||||
double volOwner = 1;
|
||||
double volNei = 1;
|
||||
double volRatio = 1;
|
||||
|
||||
fval = 1;
|
||||
|
||||
// print_patch_data(pd);
|
||||
|
||||
// std::cout << "elem_index: " << elem_index << std::endl;
|
||||
|
||||
std::vector<size_t> neighbourElements;
|
||||
pd.get_element_to_element_indices(elem_index,neighbourElements,err);
|
||||
MSQ_ERRZERO(err);
|
||||
|
||||
// std::cout << "neighbourElements: ";
|
||||
// for (size_t cellI=0; cellI < neighbourElements.size(); cellI++)
|
||||
// {
|
||||
// std::cout << neighbourElements[cellI] << " ";
|
||||
// }
|
||||
// std::cout << std::endl;
|
||||
|
||||
// std::vector<size_t> neighbourElementsNew;
|
||||
// EntityTopology entity = ownElement.get_element_type();
|
||||
// size_t connectDimension = TopologyInfo::dimension(entity)-1;
|
||||
// // 1 = connected by lines for 2D elements; 2 = connected by faces for 3D elements
|
||||
// pd.get_adjacent_entities_via_n_dim(connectDimension,elem_index,neighbourElementsNew,err);
|
||||
// MSQ_ERRZERO(err);
|
||||
//
|
||||
// std::cout << "neighbourElementsNew: ";
|
||||
// for (size_t cellI=0; cellI < neighbourElementsNew.size(); cellI++)
|
||||
// {
|
||||
// std::cout << neighbourElementsNew[cellI] << " ";
|
||||
// }
|
||||
// std::cout << std::endl;
|
||||
|
||||
// compute positive volume of owner cell
|
||||
// volOwner = ownElement.compute_unsigned_area(pd, err);
|
||||
volOwner = ownElement.compute_signed_area(pd, err);
|
||||
// MSQ_ERRRTN(err);
|
||||
|
||||
for (size_t cellI=0; cellI < neighbourElements.size(); cellI++)
|
||||
{
|
||||
MsqMeshEntity& neiElement = pd.element_by_index(neighbourElements[cellI]);
|
||||
// volNei = neiElement.compute_unsigned_area(pd, err);
|
||||
volNei = neiElement.compute_signed_area(pd, err);
|
||||
// MSQ_ERRRTN(err);
|
||||
volRatio = volOwner/(volNei+1e-15); //+SMALL
|
||||
if (fabs(volRatio) < 1.0) {volRatio = 1.0/(volRatio+1e-15);}
|
||||
fval = std::max(fval,volRatio);
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
|
@ -0,0 +1,74 @@
|
|||
/* *****************************************************************
|
||||
MESQUITE -- The Mesh Quality Improvement Toolkit
|
||||
|
||||
Copyright 2011 Oliver Borm
|
||||
|
||||
This library is free software; you can redistribute it and/or
|
||||
modify it under the terms of the GNU Lesser General Public
|
||||
License as published by the Free Software Foundation; either
|
||||
version 2.1 of the License, or (at your option) any later version.
|
||||
|
||||
This library is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
||||
Lesser General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU Lesser General Public License
|
||||
(lgpl.txt) along with this library; if not, write to the Free Software
|
||||
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
|
||||
|
||||
oli.borm@web.de
|
||||
|
||||
***************************************************************** */
|
||||
// -*- Mode : c++; tab-width: 3; c-tab-always-indent: t; indent-tabs-mode: nil; c-basic-offset: 3 -*-
|
||||
|
||||
/*! \file VolumeRatioQualityMetric.hpp
|
||||
\brief
|
||||
Header file for the Mesquite::VolumeRatioQualityMetric class
|
||||
|
||||
\author Oliver Borm
|
||||
\date 2011-01-17
|
||||
*/
|
||||
|
||||
|
||||
#ifndef VolumeRatioQualityMetric_hpp
|
||||
#define VolumeRatioQualityMetric_hpp
|
||||
|
||||
#include "Mesquite.hpp"
|
||||
#include "MsqError.hpp"
|
||||
#include "ElementQM.hpp"
|
||||
namespace MESQUITE_NS
|
||||
{
|
||||
/*! \class VolumeRatioQualityMetric
|
||||
\brief Object for computing the maximum volume ratio of
|
||||
an element, by checking all neighbour elements.
|
||||
*/
|
||||
class VolumeRatioQualityMetric : public ElementQM
|
||||
{
|
||||
public:
|
||||
VolumeRatioQualityMetric() {}
|
||||
|
||||
//! virtual destructor ensures use of polymorphism during destruction
|
||||
virtual ~VolumeRatioQualityMetric()
|
||||
{}
|
||||
|
||||
virtual std::string get_name() const;
|
||||
|
||||
int get_negate_flag() const;
|
||||
|
||||
bool evaluate( PatchData& pd,
|
||||
size_t handle,
|
||||
double& value,
|
||||
MsqError& err );
|
||||
|
||||
private:
|
||||
std::vector<Vector3D> vert;
|
||||
};
|
||||
|
||||
|
||||
} //namespace
|
||||
|
||||
|
||||
#endif // VolumeRatioQualityMetric_hpp
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,6 @@
|
|||
These are additional changes to the original mesquite src code, with two adapted Mesh Quality metrics:
|
||||
|
||||
- AspectRatioGamma is now also working for Prism, Hexahedra and Pyramids
|
||||
- VolumeRatio is a new metric, which computes the relative volume change between two adjacent cells
|
||||
|
||||
In order to use these source codes, you'll need to copy these files in a proper place of your local mesquite copy, reconfigure, recompile and install the mesquite library as well the header file in a proper place. A simple compile.sh script is shipped with this additional files, maybe it is helpfull.
|
||||
21
src/dynamicMesh/meshMotion/mesquiteMotionSolver/msqAdditonalSrc/compile.sh
Executable file
21
src/dynamicMesh/meshMotion/mesquiteMotionSolver/msqAdditonalSrc/compile.sh
Executable file
|
|
@ -0,0 +1,21 @@
|
|||
#!/bin/bash
|
||||
|
||||
rm Makefile
|
||||
rm Makefile.in
|
||||
autoreconf -fi
|
||||
|
||||
./configure \
|
||||
--enable-shared \
|
||||
--disable-imesh \
|
||||
--disable-igeom \
|
||||
--disable-irel \
|
||||
--without-cppunit \
|
||||
--enable-trap-fpe \
|
||||
--disable-function-timers \
|
||||
--enable-api-doc=HTML \
|
||||
--enable-user-guide=PDF
|
||||
|
||||
make mostlyclean-generic
|
||||
make
|
||||
|
||||
# make DESTDIR=/fu/bar install
|
||||
Reference in a new issue