Primitive patch to store face area vectors. Immersed boundary

src/foam/meshes/primitiveMesh/PrimitivePatchTemplate/PrimitivePatchClear.C

@@ -48,6 +48,7 @@ clearGeom()
 
     deleteDemandDrivenData(localPointsPtr_);
     deleteDemandDrivenData(faceCentresPtr_);
+    deleteDemandDrivenData(faceAreasPtr_);
     deleteDemandDrivenData(faceNormalsPtr_);
     deleteDemandDrivenData(pointNormalsPtr_);
 }

src/foam/meshes/primitiveMesh/PrimitivePatchTemplate/PrimitivePatchMeshData.C

@@ -376,43 +376,42 @@ template
 >
 void
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
-calcFaceNormals() const
+calcFaceAreas() const
 {
     if (debug)
     {
         Pout<< "PrimitivePatch<Face, FaceList, PointField, PointType>::"
-               "calcFaceNormals() : "
-               "calculating faceNormals in PrimitivePatch"
+               "calcFaceAreas() : "
+               "calculating faceAreas in PrimitivePatch"
             << endl;
     }
 
-    // It is considered an error to attempt to recalculate faceNormals
+    // It is considered an error to attempt to recalculate faceAreas
     // if they have already been calculated.
-    if (faceNormalsPtr_)
+    if (faceAreasPtr_)
     {
         FatalErrorIn
         (
             "PrimitivePatch<Face, FaceList, PointField, PointType>::"
-            "calcFaceNormals()"
-        )   << "faceNormalsPtr_ already allocated"
+            "calcFaceAreas()"
+        )   << "faceAreasPtr_ already allocated"
             << abort(FatalError);
     }
 
-    faceNormalsPtr_ = new Field<PointType>(this->size());
+    faceAreasPtr_ = new Field<PointType>(this->size());
 
-    Field<PointType>& n = *faceNormalsPtr_;
+    Field<PointType>& n = *faceAreasPtr_;
 
     forAll (n, faceI)
     {
         n[faceI] = this->operator[](faceI).normal(points_);
-        n[faceI] /= mag(n[faceI]) + VSMALL;
     }
 
     if (debug)
     {
         Pout<< "PrimitivePatch<Face, FaceList, PointField, PointType>::"
-               "calcFaceNormals() : "
-               "finished calculating faceNormals in PrimitivePatch"
+               "calcFaceAreas() : "
+               "finished calculating faceAreas in PrimitivePatch"
             << endl;
     }
 }

src/foam/meshes/primitiveMesh/PrimitivePatchTemplate/PrimitivePatchTemplate.C

@@ -77,6 +77,7 @@ PrimitivePatch
     localPointsPtr_(NULL),
     localPointOrderPtr_(NULL),
     faceCentresPtr_(NULL),
+    faceAreasPtr_(NULL),
     faceNormalsPtr_(NULL),
     pointNormalsPtr_(NULL)
 {
@@ -115,6 +116,7 @@ PrimitivePatch
     localPointsPtr_(NULL),
     localPointOrderPtr_(NULL),
     faceCentresPtr_(NULL),
+    faceAreasPtr_(NULL),
     faceNormalsPtr_(NULL),
     pointNormalsPtr_(NULL)
 {}
@@ -151,6 +153,7 @@ PrimitivePatch
     localPointsPtr_(NULL),
     localPointOrderPtr_(NULL),
     faceCentresPtr_(NULL),
+    faceAreasPtr_(NULL),
     faceNormalsPtr_(NULL),
     pointNormalsPtr_(NULL)
 {}
@@ -508,6 +511,26 @@ faceCentres() const
 }
 
 
+// Note: avoiding name clash
+template
+<
+    class Face,
+    template<class> class FaceList,
+    class PointField,
+    class PointType
+>
+const Foam::Field<PointType>&
+Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::areas() const
+{
+    if (!faceAreasPtr_)
+    {
+        calcFaceAreas();
+    }
+
+    return *faceAreasPtr_;
+}
+
+
 template
 <
     class Face,
@@ -521,7 +544,9 @@ faceNormals() const
 {
     if (!faceNormalsPtr_)
     {
-        calcFaceNormals();
+        // Note: stabilisation of face normals
+        faceNormalsPtr_ =
+            new Field<PointType>(areas()/(mag(areas()) + VSMALL));
     }
 
     return *faceNormalsPtr_;

src/foam/meshes/primitiveMesh/PrimitivePatchTemplate/PrimitivePatchTemplate.H

@@ -169,6 +169,9 @@ private:
         //- Face centres
         mutable Field<PointType>* faceCentresPtr_;
 
+        //- Face areas
+        mutable Field<PointType>* faceAreasPtr_;
+
         //- Face unit normals
         mutable Field<PointType>* faceNormalsPtr_;
 
@@ -214,8 +217,8 @@ private:
         //- Calculate face centres
         void calcFaceCentres() const;
 
-        //- Calculate unit face normals
-        void calcFaceNormals() const;
+        //- Calculate face areas
+        void calcFaceAreas() const;
 
         //- Calculate unit point normals
         void calcPointNormals() const;
@@ -371,6 +374,9 @@ public:
                 //- Return face centres for patch
                 const Field<PointType>& faceCentres() const;
 
+                //- Return face areas for patch.  Avoiding name clash
+                const Field<PointType>& areas() const;
+
                 //- Return face normals for patch
                 const Field<PointType>& faceNormals() const;