Merge remote-tracking branch 'origin/feature/mesquiteHexPrismSupport' into nextRelease

Conflicts:
	src/dynamicMesh/meshMotion/mesquiteMotionSolver/mesquiteMotionSolver.C
	src/dynamicMesh/meshMotion/mesquiteMotionSolver/mesquiteMotionSolver.H
This commit is contained in:
Dominik Christ 2013-07-08 13:51:02 +01:00
commit e83d2551fc
33 changed files with 20499 additions and 424 deletions

8
.gitignore vendored
View file

@ -173,4 +173,12 @@ src/lduSolvers/amg/amgPolicy/samgPolicy.H
src/lduSolvers/amg/amgPolicy/aamgPolicy.C src/lduSolvers/amg/amgPolicy/aamgPolicy.C
src/lduSolvers/amg/amgPolicy/aamgPolicy.H src/lduSolvers/amg/amgPolicy/aamgPolicy.H
# The following files are blacklisted because of a DMCA complaint by ANSYS.
src/lduSolvers/tools/PriorityArray.C
src/lduSolvers/tools/PriorityArray.H
src/lduSolvers/amg/amgPolicy/samgPolicy.C
src/lduSolvers/amg/amgPolicy/samgPolicy.H
src/lduSolvers/amg/amgPolicy/aamgPolicy.C
src/lduSolvers/amg/amgPolicy/aamgPolicy.H
# end-of-file # end-of-file

View file

@ -42,9 +42,16 @@ SourceFiles
#include "Map.H" #include "Map.H"
#include "Switch.H" #include "Switch.H"
#include "edgeList.H" #include "edgeList.H"
#include "labelPair.H"
#include "pointIOField.H" #include "pointIOField.H"
#include "MeshObject.H" #include "MeshObject.H"
#include "HashSet.H" #include "HashSet.H"
// Have gcc ignore certain warnings while including mesquite headers
#if defined(__GNUC__) && !defined(__INTEL_COMPILER)
# pragma GCC diagnostic ignored "-Wold-style-cast"
#endif
#include "Mesquite_all_headers.hpp" #include "Mesquite_all_headers.hpp"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -113,13 +120,16 @@ class mesquiteMotionSolver
scalar relax_; scalar relax_;
//- Vertex coordinate array passed into Mesquite //- Vertex coordinate array passed into Mesquite
mutable double* vtxCoords_; mutable List<double> vtxCoords_;
//- Connectivity array passed into mesquite //- Connectivity array passed into mesquite
mutable unsigned long* cellToNode_; mutable List<unsigned long> cellToNode_;
//- Flag array for vertices (fixed/free) //- Flag array for vertices (fixed/free)
mutable int* fixFlags_; mutable List<int> fixFlags_;
//- Array with element types
mutable List<Mesquite::EntityTopology> mixedTypes_;
//- Reference points //- Reference points
mutable pointIOField refPoints_; mutable pointIOField refPoints_;
@ -127,8 +137,8 @@ class mesquiteMotionSolver
//- Original points prior to surface-smoothing //- Original points prior to surface-smoothing
mutable pointField origPoints_; mutable pointField origPoints_;
//- The quality metric //- Quality metric list
word qMetric_; wordList qMetrics_;
//- Pointers to base element metric class //- Pointers to base element metric class
HashTable<autoPtr<Mesquite::QualityMetric> > qMetricTable_; HashTable<autoPtr<Mesquite::QualityMetric> > qMetricTable_;
@ -155,7 +165,9 @@ class mesquiteMotionSolver
labelList procIndices_; labelList procIndices_;
scalarField pointFractions_; scalarField pointFractions_;
labelListList procPointLabels_; labelListList procPointLabels_;
List<List<labelPair> > globalProcPoints_;
List<Map<label> > unmatchedRecvPoints_;
List<vectorField> sendSurfFields_, recvSurfFields_; List<vectorField> sendSurfFields_, recvSurfFields_;
List<Map<label> > sendSurfPointMap_, recvSurfPointMap_; List<Map<label> > sendSurfPointMap_, recvSurfPointMap_;
@ -207,7 +219,7 @@ class mesquiteMotionSolver
// Component-wise sumMag // Component-wise sumMag
scalar cmptSumMag(const vectorField& field); scalar cmptSumMag(const vectorField& field);
// Transfer buffers after divergence compute. // Transfer buffers for surface point fields
void transferBuffers(vectorField &field); void transferBuffers(vectorField &field);
// Apply boundary conditions // Apply boundary conditions
@ -283,9 +295,6 @@ class mesquiteMotionSolver
scalar& beta scalar& beta
); );
// Clear out addressing
void clearOut();
// Update on topology change // Update on topology change
void update(const mapPolyMesh&); void update(const mapPolyMesh&);
@ -341,4 +350,3 @@ public:
#endif #endif
// ************************************************************************* // // ************************************************************************* //

View file

@ -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 &centroid, 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

View file

@ -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;
}

View file

@ -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

View file

@ -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;
}

View file

@ -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

View file

@ -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.

View 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

View file

@ -0,0 +1,42 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.5-dev |
| \\ / A nd | Revision: 1557 |
| \\/ M anipulation | Web: http://www.OpenFOAM.org |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
location "0";
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
referenceLevel (0 0 0);
boundaryField
{
bottomWall
{
type fixedValue;
value uniform (0 0 0);
}
topWall
{
type fixedValue;
value uniform (0 0 0);
}
sideWall
{
type fixedValue;
value uniform (0 0 0);
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View file

@ -0,0 +1,61 @@
/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.4.1-dev |
| \\ / A nd | Web: http://www.openfoam.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
root "";
case "";
instance "";
local "";
class tetPointVectorField;
object motionU;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
bottomWall
{
type fixedValue;
value uniform (0 0 0);
}
sideWall
{
type zeroGradient; //slip;
}
topWall
{
type oscillatingFixedValue;
refValue uniform (0 0 0);
amplitude uniform (0 0 0.5);
frequency 0.1;
//type fixedValue;
//value uniform (0 0 0.1);
}
outlet
{
type fixedValue;
value uniform (0 0 0);
}
}
// ************************************************************************* //

View file

@ -0,0 +1,44 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.5-dev |
| \\ / A nd | Revision: 1557 |
| \\/ M anipulation | Web: http://www.OpenFOAM.org |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
referenceLevel 0;
boundaryField
{
bottomWall
{
type fixedValue;
value uniform 0.0;
}
topWall
{
type fixedValue;
value uniform 0.0;
}
sideWall
{
type fixedValue;
value uniform 0.0;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View file

@ -0,0 +1,6 @@
#!/bin/sh
. $WM_PROJECT_DIR/bin/tools/CleanFunctions
cleanTimeDirectories
rm -rf processor* > /dev/null 2>&1

View file

@ -0,0 +1,10 @@
#!/bin/bash
# Source tutorial run functions
. $WM_PROJECT_DIR/bin/tools/RunFunctions
# Get application name
application="moveDynamicMesh"
runApplication decomposePar
runParallel $application 4

View file

@ -0,0 +1,122 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.5-dev |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object dynamicMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Select the type of dynamicFvMesh, and load custom library
dynamicFvMesh dynamicMotionSolverFvMesh;
solver mesquiteMotionSolver;
//solver laplaceFaceDecomposition;
//diffusivity uniform 1.0;
//diffusivity quadratic;
//distancePatches (topWall);
//frozenDiffusion off;
mesquiteOptions
{
// Optimization metric
// optMetric AspectRatioGamma;
optMetric
{
// firstMetric EdgeLength;
firstMetric MeanRatio;
// secondMetric MeanRatio;
secondMetric AspectRatioGamma;
}
// Objective function
objFunction CompositeOFAdd;
// Optimization algorithm
optAlgorithm FeasibleNewton;
// optAlgorithm ConjugateGradient;
// Termination criteria sub-dictionary (takes default values if not specified)
// Specifying an empty sub-dictionary terminates with available options
tcInner
{
//relGradL2 1e-2;
//cpuTime 0.5;
iterationLimit 5;
}
// tcOuter
// {}
sliverThreshold 0.05;
// For composite functions, two objectives need to be specified
firstFunction LInf;
secondFunction LPtoP;
// For scaled functions, scale and objective needs to be specified
// scaleFunction PMeanP;
// scale 1.5;
// Power value for the LPtoP objective function
pValue 20;
power 2;
// Specify a tolerance for the CG solver
tolerance 1e-2;
// Specify number of CG sweeps
nSweeps 1;
// Specify a relaxation factor, if necessary
//relaxationFactor 0.1;
// Specify slip patches for the motionSolver
slipPatches
{
sideWall;
topWall;
bottomWall;
}
cylindricalConstraints
{
// Specify options per slip patch
sideWall
{
axisPoint (0.0 0.0 0.0);
axisVector (0.0 0.0 1.0);
radius 1.0;
}
}
// Specify fixedValue patches for the motionSolver
fixedValuePatches
{
topWall
{
//type angularOscillatingDisplacement;
//amplitude -0.0125;
type oscillatingDisplacement;
amplitude (0 0 0.01);
axis (1 0 0);
origin (0 0 3);
angle0 0.0;
omega 0.15;
value uniform (0 0 0);
}
}
// Specify interval for surface smoothing
surfInterval 1;
}
// ************************************************************************* //

View file

@ -0,0 +1,40 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM Extend Project: Open source CFD |
| \\ / O peration | Version: 1.6-ext |
| \\ / A nd | Web: www.extend-project.de |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class polyBoundaryMesh;
location "constant/polyMesh";
object boundary;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
3
(
bottomWall
{
type wall;
nFaces 89;
startFace 3676;
}
topWall
{
type wall;
nFaces 89;
startFace 3765;
}
sideWall
{
type wall;
nFaces 480;
startFace 3854;
}
)
// ************************************************************************* //

View file

@ -0,0 +1,21 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM Extend Project: Open source CFD |
| \\ / O peration | Version: 1.6-ext |
| \\ / A nd | Web: www.extend-project.de |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class regIOobject;
location "constant/polyMesh";
object cellZones;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
0
()
// ************************************************************************* //

View file

@ -0,0 +1,21 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM Extend Project: Open source CFD |
| \\ / O peration | Version: 1.6-ext |
| \\ / A nd | Web: www.extend-project.de |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class regIOobject;
location "constant/polyMesh";
object faceZones;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
0
()
// ************************************************************************* //

File diff suppressed because it is too large Load diff

File diff suppressed because it is too large Load diff

File diff suppressed because it is too large Load diff

View file

@ -0,0 +1,21 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM Extend Project: Open source CFD |
| \\ / O peration | Version: 1.6-ext |
| \\ / A nd | Web: www.extend-project.de |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class regIOobject;
location "constant/polyMesh";
object pointZones;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
0
()
// ************************************************************************* //

File diff suppressed because it is too large Load diff

View file

@ -0,0 +1,32 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM Extend Project: Open source CFD |
| \\ / O peration | Version: 1.6-ext |
| \\ / A nd | Web: www.extend-project.de |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class wordList;
location "constant/polyMesh";
object zoneToPatchName;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
8
(
unknown
unknown
unknown
bottomWall
topWall
sideWall
unknown
default-interior
)
// ************************************************************************* //

View file

@ -0,0 +1,30 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.4 |
| \\ / A nd | Web: http://www.openfoam.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
root "";
case "";
instance "";
local "";
class dictionary;
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu nu [0 2 -1 0 0 0 0] 0.01;
transportModel Newtonian;
// ************************************************************************* //

View file

@ -0,0 +1,62 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.5-dev |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
startTime 0;
startFrom startTime;
endTime 50.0;
stopAt endTime;
deltaT 0.1;
adjustTimeStep no;
maxCo 0.8;
writeControl timeStep;
writeInterval 1;
writeFormat ascii;
writeCompression uncompressed;
timeFormat general;
//timePrecision 15;
//writePrecision 11;
runTimeModifiable yes;
//functions
//(
// meshFluxes
// {
// Type of functionObject
// type meshFluxes;
// Where to load it from (if not already in solver)
// functionObjectLibs ("libmeshFluxes.so");
// }
//);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View file

@ -0,0 +1,46 @@
/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.3 |
| \\ / A nd | Web: http://www.openfoam.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
root "";
case "";
instance "";
local "";
class dictionary;
object decomposeParDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
numberOfSubdomains 4;
//method simple;
//method metis;
method parMetis;
tolerance 1e-3;
distributed no;
simpleCoeffs
{
n (1 1 1);
delta 0.001;
}
roots
(
);
// ************************************************************************* //

View file

@ -0,0 +1,70 @@
/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.3 |
| \\ / A nd | Web: http://www.openfoam.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
root "";
case "";
instance "";
local "";
class dictionary;
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
ddtSchemes
{
default Euler;
}
gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss linear;
}
laplacianSchemes
{
default none;
laplacian(nu,U) Gauss linear corrected;
laplacian(rUA,pcorr) Gauss linear corrected;
laplacian(rUA,p) Gauss linear corrected;
}
interpolationSchemes
{
default linear;
interpolate(HbyA) linear;
interpolate(1|A) linear;
}
snGradSchemes
{
default corrected;
}
fluxRequired
{
default no;
pcorr;
p;
}
// ************************************************************************* //

View file

@ -0,0 +1,91 @@
/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.3 |
| \\ / A nd | Web: http://www.openfoam.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
root "";
case "";
instance "";
local "";
class dictionary;
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers
{
pcorr
{
solver PCG;
preconditioner
{
type DIC;
}
minIter 0;
maxIter 1000;
tolerance 1e-02;
relTol 0;
};
p
{
solver PCG;
preconditioner
{
type DIC;
}
minIter 0;
maxIter 1000;
tolerance 1e-06;
relTol 0.05;
};
pFinal
{
solver PCG;
preconditioner
{
type DIC;
}
minIter 0;
maxIter 1000;
tolerance 1e-06;
relTol 0;
};
U
{
solver PCG;
preconditioner
{
type DILU;
}
minIter 0;
maxIter 1000;
tolerance 1e-05;
relTol 0;
};
}
PISO
{
nCorrectors 2;
nNonOrthogonalCorrectors 0;
}
// ************************************************************************* //

View file

@ -0,0 +1,30 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.5-dev |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object tetFemSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers
{
motionU
{
solver GMRES;
nDirections 8;
preconditioner Cholesky;
tolerance 1e-09;
relTol 0.0001;
}
}
// ************************************************************************* //