slim3d Forcing Zero at Bnd now availabl on slim3d

It is forced on both 2D and 3D variables

modules/slim3d/slim3dSolver.cpp

@@ -90,7 +90,11 @@ slim3dSolver::slim3dSolver(const std::string meshFile3d, std::vector<std::string
 
   copy2d3d = new dgSW3dDofCopy(columnInfo);
   verticalBottomRemover = new dgSW3dVerticalBottomRemover(columnInfo);
-  forceZeroUVAtBnd = new slim3dForceZeroAtBnd(groups3d, physicalTags2d);
+  std::vector<std::string> zeroBndTags;
+  zeroBndTags.push_back("coast");
+  zeroBndTags.push_back("vertical_bottom");
+  forceZeroUVAtBnd = new slim3dForceZeroAtBnd(groups3d, zeroBndTags);
+  forceZeroUVAtBnd2d = new slim3dForceZeroAtBnd(groups2d, zeroBndTags);
 }
 
 slim3dSolver::~slim3dSolver()
@@ -103,6 +107,7 @@ slim3dSolver::~slim3dSolver()
   delete _equations;
   delete copy2d3d;
   delete forceZeroUVAtBnd;
+  delete forceZeroUVAtBnd2d;
   if ( wSolver ) delete wSolver;
   //if ( depthIntSolver ) delete depthIntSolver;
   if ( rSolver ) delete rSolver;

modules/slim3d/slim3dSolver.h

@@ -78,6 +78,7 @@ public:
   dgCGMeanFilter *dgCG3d, *dgCG2d;
 //   dgCGStructure *dgCG3d, *dgCG2d;
   slim3dForceZeroAtBnd *forceZeroUVAtBnd;
+  slim3dForceZeroAtBnd *forceZeroUVAtBnd2d;
 
   /** Create a solver instance and load 3d and 2d mesh */
   slim3dSolver(const std::string meshFile3d, std::vector<std::string> bottomBndTags = std::vector<std::string>(), std::vector<std::string> topBndTags = std::vector<std::string>(), std::string periodicMapFile = "");

modules/slim3d/slim3dTimeIntegrator.cpp

@@ -1128,6 +1128,8 @@ void slim3dTimeIntegratorPC::advanceOneTimeStep()
 //        uvAvDof2d_old->copy(*d->uvAvDof2d);
         e->uv2dEq->computeAllTerms(time2d,*d->uvAvDof2d,*uvAvDof2d_K,false);
         uvAvDof2d_K->multiplyByInvMassMatrix();
+        if ( s->getForceZeroUVAtBnd() )
+          s->forceZeroUVAtBnd2d->forceZeroAtBnd(uvAvDof2d_K);
         if (iterAB == 0) {
           // Euler step to begin with !
           d->uvAvDof2d->axpy(*uvAvDof2d_K,_dt);

modules/slim3d/slim3dUtils.cpp

@@ -1185,25 +1185,34 @@ void slim3dForceZeroAtBnd::forceZeroAtBnd(dgDofContainer* dof)
 {
   dgFullMatrix<double> data;
   fullMatrix<double> normals;
+  functorCache cacheMap(functorCache::NODE_MODE, _groups);
   for (int iFaceGroup = 0; iFaceGroup < _groups->getNbFaceGroups(); ++iFaceGroup){
     dgGroupOfFaces *faceGroup= _groups->getFaceGroup(iFaceGroup);
-    std::vector<std::string>::iterator it = std::find(_physicalTags.begin(), _physicalTags.end(), faceGroup->physicalTag()); 
+    std::vector<std::string>::iterator it = std::find(_physicalTags.begin(), _physicalTags.end(), faceGroup->physicalTag());
     if ( it != _physicalTags.end() ){
       for (size_t iFace = 0; iFace < faceGroup->size(); ++iFace){
         const std::vector<int> &cl = faceGroup->closure(iFace, 0);
         size_t iGroup = _groups->getElementGroupId(&faceGroup->elementGroup(0));
         size_t iElem = faceGroup->elementId(iFace, 0);
         dof->getElementProxy(iGroup, iElem, data);
-        faceGroup->normal(*_jacobians, iFace, 0, normals);
+        //printf("%s, %d ..\n", faceGroup->physicalTag().c_str(), iFace);
+        //faceGroup->normal(*_jacobians, iFace, 0, normals);
+        const dgGroupOfElements *groupIn = &faceGroup->elementGroup(0);
+        cacheMap.setGroup(groupIn);
+        const dgMeshJacobian & jacobians = cacheMap.jacobian();
+        faceGroup->normal(jacobians, iFace, 0, normals);
         for (int iNode = 0; iNode < faceGroup->getNbNodes(); ++iNode){
           if(_alreadyForced[iGroup][iElem][cl[iNode]]){
             data(cl[iNode],0) = 0;
             data(cl[iNode],1) = 0;
           }
           else {
-            double dataN = data(cl[iNode], 0) * normals(iNode, 0) + data(cl[iNode], 1) * normals(iNode, 1);
-            data(cl[iNode], 0) = data(cl[iNode], 0) - normals(iNode, 0) * dataN;  
-            data(cl[iNode], 1) = data(cl[iNode], 1) - normals(iNode, 1) * dataN;  
+            //double dataN = data(cl[iNode], 0) * normals(iNode, 0) + data(cl[iNode], 1) * normals(iNode, 1);
+            //data(cl[iNode], 0) = data(cl[iNode], 0) - normals(iNode, 0) * dataN;  
+            //data(cl[iNode], 1) = data(cl[iNode], 1) - normals(iNode, 1) * dataN;  
+            double dataN = data(cl[iNode], 0) * normals(0, 0) + data(cl[iNode], 1) * normals(0, 1);
+            data(cl[iNode], 0) = data(cl[iNode], 0) - normals(0, 0) * dataN;  
+            data(cl[iNode], 1) = data(cl[iNode], 1) - normals(0, 1) * dataN;  
           }
         }
       }