Commit 69e8d27b authored by Toni Mikael Lassila's avatar Toni Mikael Lassila Committed by Antonio Cervone
Browse files

Added method for computing approximate normal vector on a boundary surface

parent 42d9914d
......@@ -183,6 +183,19 @@ public:
template< typename VectorType >
Vector average( const VectorType& field, const markerID_Type& flag, UInt feSpace = 0, UInt nDim = 1);
This method computes an average of a scalar field on the boundary section "flag"
@ingroup boundary_methods
\tparam VectorType Vector type. Basic policy for type VectorType: operator[] available
\param field is intended to be a vector or a scalar (pressure in NS problem),
this method computes the average value of field on section "flag"
\return the approximate normal vector
Vector normal( const markerID_Type& flag, UInt feSpace = 0, UInt nDim = nDimensions);
#if 0
......@@ -789,6 +802,70 @@ Vector PostProcessingBoundary<MeshType>::average( const VectorType& field, const
return fieldAverage / measure;
// approximate normal for a certain marker
template<typename MeshType>
Vector PostProcessingBoundary<MeshType>::normal( const markerID_Type& flag, UInt feSpace, UInt nDim )
// Each processor computes the normal on his own flagged facets --> normalScatter
// At the end I'll reduce the process normals --> normal
Vector normalScatter(3), normal(3);
// I need the global DOF ID to query the vector
// dofVectorIndex is the index of the dof in the data structure of PostProcessingBoundary class
// dofGlobalId is the corresponding ID in the GLOBAL mesh (prior to partitioning)
UInt dofVectorIndex, dofGlobalId;
// list of flagged facets on current processor
std::list<ID> facetList( M_boundaryMarkerToFacetIdMap[flag] );
typedef std::list<ID>::iterator Iterator;
// Nodal values of field in the current facet
Vector localFieldVector(nDim * M_numTotalDofPerFacetVector[feSpace]);
// Loop on facetList
for (Iterator j=facetList.begin(); j != facetList.end(); ++j)
// Updating quadrature data on the current facet
M_currentBdFEPtrVector[feSpace]->updateMeasNormalQuadPt(M_meshPtr->boundaryFacet(*j) );
// Quadrature formula
// Loop on quadrature points
for (UInt iq=0; iq< M_currentBdFEPtrVector[feSpace]->nbQuadPt(); ++iq)
// Dot product
// Loop on components
for (UInt iComponent =0; iComponent<nDim; ++iComponent)
// Interpolation
// Loop on local dof
for (ID iDof=0; iDof<M_numTotalDofPerFacetVector[feSpace]; ++iDof)
// Extracting nodal values of field in the current facet
dofVectorIndex = M_vectorNumberingPerFacetVector[feSpace][ ( UInt ) *j ][ iDof ];
dofGlobalId = M_dofGlobalIdVector[feSpace][dofVectorIndex]; // this is in the GLOBAL mesh
normalScatter(iComponent) += M_currentBdFEPtrVector[feSpace]->weightMeas(iq)
* M_currentBdFEPtrVector[feSpace]->phi(Int(iDof),iq)
* M_currentBdFEPtrVector[feSpace]->normal(Int(iComponent),iq);
// Reducing per-processor values
M_epetraMapPtr->comm().SumAll( &normalScatter(0), &normal(0), 1 );
M_epetraMapPtr->comm().SumAll( &normalScatter(1), &normal(1), 1 );
M_epetraMapPtr->comm().SumAll( &normalScatter(2), &normal(2), 1 );
// Scale normal to unity length
Real nn = sqrt(pow(normal(0),2) + pow(normal(1),2) + pow(normal(2),2));
return (normal / nn);
// Measure of patches on the boundary
template<typename MeshType>
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment