Merge branch 'fix_tangent_bbnd' into 'master'

fix tangential vector for BBND elements See merge request !215
NGSolve · Aug 9, 2017 · 0703290 · 0703290
2 parents df5493b + 1fb555c
commit 0703290
Show file tree

Hide file tree

Showing 2 changed files with 9 additions and 44 deletions.
diff --git a/fem/intrule.cpp b/fem/intrule.cpp
@@ -136,48 +136,6 @@ namespace ngfem
     return eltrans->VB();
   }
 
-
-  template <int S, int R, typename SCAL>
-  MappedIntegrationPoint<S,R,SCAL> :: 
-  MappedIntegrationPoint (const IntegrationPoint & aip,
-			  const ElementTransformation & aeltrans)
-    : DimMappedIntegrationPoint<R,SCAL> (aip, aeltrans)
-  {
-    this->eltrans->CalcPointJacobian(this->IP(), this->point, dxdxi);
-
-    if (S == R)
-      {
-	det = Det (dxdxi); 
-	if(det == 0) 
-	  { 
-	    throw Exception ("Jacobi-det is zero");
-	  } 
-	if (det < 0 && 0)
-	  {
-	    throw Exception ("Jacobi-det is negative");
-	  }
-	// dxidx = Inv (dxdxi);
-      }
-    else
-      {
-	if (R == 3)
-	  {
-            normalvec = Cross (Vec<3,SCAL> (dxdxi.Col(0)),
-                               Vec<3,SCAL> (dxdxi.Col(1)));
-            det = L2Norm (normalvec);
-            normalvec /= det;
-	  }
-	else
-	  {
-	    det = sqrt ( sqr (dxdxi(0,0)) + sqr (dxdxi(1,0)));
-
-	    normalvec(0) = -dxdxi(1,0) / det;
-	    normalvec(1) = dxdxi(0,0) / det;
-	  }
-      }
-    this->measure = fabs (det);
-  }
-
 
   template <int S, int R, typename SCAL>
   void MappedIntegrationPoint<S,R,SCAL> :: 

diff --git a/fem/intrule.hpp b/fem/intrule.hpp
@@ -368,7 +368,12 @@ namespace ngfem
     NGS_DLL_HEADER MappedIntegrationPoint () = default;
     ///
     NGS_DLL_HEADER MappedIntegrationPoint (const IntegrationPoint & aip,
-					   const ElementTransformation & aeltrans);
+					   const ElementTransformation & aeltrans)
+      : DimMappedIntegrationPoint<DIMR,SCAL> (aip, aeltrans)
+    {
+      this->eltrans->CalcPointJacobian(this->IP(), this->point, dxdxi);
+      this->Compute();
+    }
 
     INLINE MappedIntegrationPoint (const IntegrationPoint & aip,
 			    const ElementTransformation & aeltrans,
@@ -406,6 +411,7 @@ namespace ngfem
 				     Vec<3,SCAL> (dxdxi.Col(1)));
 		  det = L2Norm (normalvec);
 		  normalvec /= det;
+                  tangentialvec = TSCAL(0.0);
 		}
 	      else
 		{
@@ -422,13 +428,14 @@ namespace ngfem
 
 	      normalvec(0) = -dxdxi(1,0) / det;
 	      normalvec(1) = dxdxi(0,0) / det;
+              tangentialvec = TSCAL(0.0);
 	    }
 	  else
 	    {
 	      det = 1.0;
 	      normalvec = 1.0;
+              tangentialvec = TSCAL(0.0);
 	    }
-	  tangentialvec = TSCAL(0.0);
 	}
       this->measure = fabs (det);
     }