// ================================================================================
/*!
* \brief Creates a dual mesh from the given one
- *
- * The dual mesh consists of dual cells corresponding to nodes of the input mesh.
- * The dual cell is a holygon in 2D space and a polyhedron in 3D space. The dual
- * cell is bound by edges (2D) or triangles (3D) connecting barycentres of cells
- * around the input node with middles of edges ending at the input node.
*/
DualMESH::DualMESH(const MESH & mesh ):MESH()
{
// Make dual nodal connectivity.
// nodes of dual mesh are:
// 1. nodes at bary centres of input cells
- // 2. nodes at middles of input edges
+ // 2. nodes at middles of input edges (hereafter called "middle nodes")
// 3. nodes at bary centres of side faces of tetrahedrons (3D only)
- // 4. boundary nodes
+ // 4. nodes coincident with input boundary nodes
- _connectivity = new CONNECTIVITY( /*numberOfTypes = */0 );
+ _connectivity = new CONNECTIVITY( /*numberOfClassicTypes = */0 );
_connectivity->setEntityDimension( _spaceDimension );
list< pair< int, int > > bndNodes; // pairs of input and dual mesh boundary nodes
int nbBndNodes = 0; // counter of bndNodes
- set< TLink > edges; // links between nodes storing id of middle node
- set< TTriaFace > triangles; // triangles storing bary centre node id (3D only)
- typedef set<TLink>::iterator TEdgeIt;
+ set< TLink > edges; // input edges storing id of middle node
+ set< TTriaFace > triangles; // input triangles storing bary centre node id (3D only)
+ typedef set<TLink>::iterator TEdgeIt;
typedef set<TTriaFace>::iterator TFaceIt;
// --------------------------------------------------------------------------------
for ( iNode = 1; iNode <= nbDualCells; ++iNode )
{
// to make connectivity of poly around an input node, use map of middle node
- // of edge to TLink between the middle node and a neighbor middle node.
- // Each TLink represents two edges of a polygon: middle1-bary and bary-middle2
+ // to TLink between the middle node and a neighbor middle node.
+ // Each TLink between middle nodes represents two edges of a polygon:
+ // middle1-bary and bary-middle2
map< int, TLink > middle2BaryLink;
typedef map< int, TLink >::iterator TMid2LinkIt;
// loop on triangles around iNode
else {
node[0] = triaConn[0]; node[1] = triaConn[1];
}
- // get middle nodes on edges iNode-node[0] and iNode-node[1]
+ // get middle nodes on two edges: 1) iNode-node[0] and 2) iNode-node[1]
int middle[2];
for ( int i = 0; i < 2; ++i ) {
int nextMiddle = nbCellBary + edges.size() + nbBndNodes + 1; // next free node id
{
bc[3]=0; // for no solution
- // Find bc by solving system of 3 equations
- // Equation for X:
+ // Find bc by solving system of 3 equations using Gaussian elimination algorithm
// bc1*( x1 - x4 ) + bc2*( x2 - x4 ) + bc3*( x3 - x4 ) = px - x4
+ // bc1*( y1 - y4 ) + bc2*( y2 - y4 ) + bc3*( y3 - y4 ) = px - y4
+ // bc1*( z1 - z4 ) + bc2*( z2 - z4 ) + bc3*( z3 - z4 ) = px - z4
const int nbCol=4, nbRow=3;
double T[nbRow][nbCol]=
{ n[0][_Y]-n[3][_Y], n[1][_Y]-n[3][_Y], n[2][_Y]-n[3][_Y], p[_Y]-n[3][_Y] },
{ n[0][_Z]-n[3][_Z], n[1][_Z]-n[3][_Z], n[2][_Z]-n[3][_Z], p[_Z]-n[3][_Z] }};
- // make upper triangular matrix
+ // make upper triangular matrix (forward elimination)
int iR[nbRow] = { 0, 1, 2 };
return; // no solution
}
tRow[ 3 ] /= tRow[ 2 ];
- // calculate solution: backward substitution
+
+ // calculate solution (back substitution)
bc[ 2 ] = tRow[ 3 ];
else return false;
}
else
- if(i2next == istart2)
- return false;
- else
- {
- if(distance2<T,dim>(L1,i1next*dim, L2,i2next*dim) > epsilon )
- return false;
- else
- {
- i1 = i1next;
- i2 = i2next;
- i1next = ( i1 + 1 ) % size1;
- i2next = ( i2 + sign + size2 ) % size2;
- }
- }
+ if(i2next == istart2)
+ return false;
+ else
+ {
+ if(distance2<T,dim>(L1,i1next*dim, L2,i2next*dim) > epsilon )
+ return false;
+ else
+ {
+ i1 = i1next;
+ i2 = i2next;
+ i1next = ( i1 + 1 ) % size1;
+ i2next = ( i2 + sign + size2 ) % size2;
+ }
+ }
}
}
- /*! Tests if two list of nodes (not necessarily distincts) describe the same polygon.*/
- /*! Existence of multiple points in the list is considered.*/
- template<class T, int dim>
- bool checkEqualPolygons(T * L1, T * L2, double epsilon)
- {
- if(L1==NULL || L2==NULL)
- {
- std::cout << "Warning InterpolationUtils.hxx:checkEqualPolygonsPointer: Null pointer " << std::endl;
- throw(Exception("big error: not closed polygon..."));
- }
-
- int size1 = (*L1).size()/dim;
- int size2 = (*L2).size()/dim;
- int istart1 = 0;
- int istart2 = 0;
-
- while( istart2 < size2 && distance2<T,dim>(L1,istart1*dim, L2,istart2*dim) > epsilon ) istart2++;
-
- if(istart2 == size2)
- {
- return (size1 == 0) && (size2 == 0);
- }
- else
- return checkEqualPolygonsOneDirection<T,dim>( L1, L2, size1, size2, istart1, istart2, epsilon, 1)
- || checkEqualPolygonsOneDirection<T,dim>( L1, L2, size1, size2, istart1, istart2, epsilon, -1);
+ /*! Tests if two list of nodes (not necessarily distincts) describe the same polygon.*/
+ /*! Existence of multiple points in the list is considered.*/
+ template<class T, int dim>
+ bool checkEqualPolygons(T * L1, T * L2, double epsilon)
+ {
+ if(L1==NULL || L2==NULL)
+ {
+ std::cout << "Warning InterpolationUtils.hxx:checkEqualPolygonsPointer: Null pointer " << std::endl;
+ throw(Exception("big error: not closed polygon..."));
+ }
- }
+ int size1 = (*L1).size()/dim;
+ int size2 = (*L2).size()/dim;
+ int istart1 = 0;
+ int istart2 = 0;
+
+ while( istart2 < size2 && distance2<T,dim>(L1,istart1*dim, L2,istart2*dim) > epsilon ) istart2++;
+
+ if(istart2 == size2)
+ {
+ return (size1 == 0) && (size2 == 0);
+ }
+ else
+ return checkEqualPolygonsOneDirection<T,dim>( L1, L2, size1, size2, istart1, istart2, epsilon, 1)
+ || checkEqualPolygonsOneDirection<T,dim>( L1, L2, size1, size2, istart1, istart2, epsilon, -1);
+
+ }
}
