add P1nc (Crouziex -Raviart FE?) on meshS

Author: frederichecht

plugin/seq/Element_P1nc_3d.cpp

@@ -42,7 +42,75 @@
 /// ---------------------------------------------------------------
 namespace Fem2D {
 
-  class TypeOfFE_P1nc_3d : public TypeOfFE_Lagrange< Mesh3 > {
+
+class TypeOfFE_P1nc_S : public TypeOfFE_Lagrange< MeshS > {
+ public:
+ typedef MeshS Mesh;
+ typedef MeshS::Element Element;
+ typedef GFElement< MeshS > FElement;
+    static int dfon[];
+    static const int d = Mesh::Rd::d;
+     TypeOfFE_P1nc_S(): TypeOfFE_Lagrange<MeshS>(-1,2) {  }
+
+     void FB(const What_d whatd, const MeshS &Th, const MeshS::Element &K, const MeshS::RdHat &PHat,
+            RNMK_ &val) const;
+   };
+
+
+
+
+
+
+void TypeOfFE_P1nc_S::FB(const What_d whatd, const MeshS &Th, const MeshS::Element &K,
+                         const MeshS::RdHat &PHat, RNMK_ &val) const {
+  
+  assert(val.N( ) >= 3);    // 4 degrees of freedom
+  assert(val.M( ) == 1);     // 1 components
+/*
+ Dof  Numbering  1 dof / face   the dof a associated to a vertex face
+and the vertex a numbering in increase way.
+
+we have to numbering    î on ref element with
+
+ So we need of a permutation p to insure the compatibilit between adjacent element
+ 
+*/
+ 
+
+
+  double l[3];
+  PHat.toBary(l);
+  val=0; ;
+  RN_ f0(val('.', 0, op_id));
+  if (whatd & Fop_D0) {
+      RN_ f0(val('.',0,0));
+      f0[0] = 1-l[0]*2;
+      f0[1] = 1-l[1]*2;
+      f0[2] = 1-l[2]*2;
+ 
+
+     
+  }
+  if (whatd & (Fop_D1 )) {
+      const unsigned int iDx =whatd & Fop_dx , iDy =whatd & Fop_dy , iDz =whatd & Fop_dz;
+      R3 Dl[4];
+      K.Gradlambda(Dl);
+       int k=0;
+      RN_ f0x(val('.', 0, op_dx));
+      RN_ f0y(val('.', 0, op_dy));
+      RN_ f0z(val('.', 0, op_dz));
+      for(int i=0; i<3;++i)
+      {
+          if(iDx) f0x[i] = -Dl[i].x*2.;
+          if(iDy) f0y[i] = -Dl[i].y*2.;
+          if(iDz) f0z[i] = -Dl[i].z*2.;
+      }
+      
+  }
+
+}
+
+class TypeOfFE_P1nc_3d : public TypeOfFE_Lagrange< Mesh3 > {
    public:
    typedef Mesh3 Mesh;
    typedef Mesh3::Element Element;
@@ -115,10 +183,13 @@ void TypeOfFE_P1nc_3d::FB(const What_d whatd, const Mesh &Th, const Mesh3::Eleme
 
 }
  static TypeOfFE_P1nc_3d P1nc_3d;
+ static TypeOfFE_P1nc_S P1nc_S;
  GTypeOfFE< Mesh3 > &Elm_P1nc_3d(P1nc_3d);
+ GTypeOfFE< MeshS > &Elm_P1nc_S(P1nc_S);
+
+ static AddNewFE3 TFE_P1pnc_3d("P1nc3d", &Elm_P1nc_3d,"P1nc");
+ static AddNewFES TFE_P1pnc_S("P1ncS", &Elm_P1nc_S,"P1nc");
 
- static AddNewFE3 TFE_P2pnc_3d("P1nc3d", &Elm_P1nc_3d);
-  
 }    // namespace Fem2D
 
 // --- fin --