adding skew-symmetric integrator

fem/nonlininteg.cpp

@@ -952,4 +952,118 @@ void SkewSymmetricVectorConvectionNLFIntegrator::AssembleElementGrad(
    }
 }
 
+void ConvectiveNSENLFIntegrator::AssembleElementGrad(
+   const FiniteElement &el,
+   ElementTransformation &trans,
+   const Vector &elfun,
+   DenseMatrix &elmat)
+{
+   const int nd = el.GetDof();
+   const int dim = el.GetDim();
+
+   shape.SetSize(nd);
+   dshape.SetSize(nd, dim);
+   dshapex.SetSize(nd, dim);
+   elmat.SetSize(nd * dim);
+   pelmat_comp.SetSize(nd);
+   elmat_comp.SetSize(nd);
+   gradEF.SetSize(dim);
+
+   EF.UseExternalData(elfun.GetData(), nd, dim);
+
+   Vector vec1(dim), vec2(dim), vec3(nd);
+
+   const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, trans);
+
+   elmat = 0.0;
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+      trans.SetIntPoint(&ip);
+
+      el.CalcShape(ip, shape);
+      el.CalcDShape(ip, dshape);
+
+      double w = trans.Weight();
+      w = Q ? Q->Eval(trans, ip) * ip.weight / w : ip.weight / w;
+      const double w_non = ip.weight;
+
+      EF.MultTranspose(shape, vec1); // u^n
+
+      trans.AdjugateJacobian().Mult(vec1, vec2);
+      Mult(dshape, trans.AdjugateJacobian(), dshapex);
+
+      vec2 *= w_non;
+      dshape.Mult(vec2, vec3); // (u^n \cdot grad u^{n+1})
+      MultVWt(shape, vec3, elmat_comp); // (u^n \cdot grad u^{n+1},v)
+      Mult_a_AAt(w, dshapex, pelmat_comp); // (grad u^{n+1},v)
+
+      for (int i = 0; i < dim; i++)
+      {
+         elmat.AddMatrix(elmat_comp, i * nd, i * nd);
+         elmat.AddMatrix(pelmat_comp, i * nd, i * nd);
+      }
+   }
+}
+
+
+
+void SkewSymmetricNSENLFIntegrator::AssembleElementGrad(
+   const FiniteElement &el,
+   ElementTransformation &trans,
+   const Vector &elfun,
+   DenseMatrix &elmat)
+{
+   const int nd = el.GetDof();
+   const int dim = el.GetDim();
+
+   shape.SetSize(nd);
+   dshape.SetSize(nd, dim);
+   dshapex.SetSize(nd, dim);
+   elmat.SetSize(nd * dim);
+   pelmat_comp.SetSize(nd);
+   elmat_comp.SetSize(nd);
+   gradEF.SetSize(dim);
+
+   DenseMatrix elmat_comp_T(nd);
+
+   EF.UseExternalData(elfun.GetData(), nd, dim);
+
+   Vector vec1(dim), vec2(dim), vec3(nd), vec4(dim), vec5(nd);
+
+   const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, trans);
+
+   elmat = 0.0;
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+      trans.SetIntPoint(&ip);
+
+      el.CalcShape(ip, shape);
+      el.CalcDShape(ip, dshape);
+
+      double w = trans.Weight();
+      w = Q ? Q->Eval(trans, ip) * ip.weight / w : ip.weight / w;
+      const double w_non = ip.weight;
+
+      EF.MultTranspose(shape, vec1); // u^n
+
+      trans.AdjugateJacobian().Mult(vec1, vec2);
+      Mult(dshape, trans.AdjugateJacobian(), dshapex);
+
+      vec2 *= w_non;
+      dshape.Mult(vec2, vec3); // (u^n \cdot grad u^{n+1})
+      MultVWt(shape, vec3, elmat_comp); // (u^n \cdot grad u^{n+1},v)
+      elmat_comp_T.Transpose(elmat_comp);
+      Mult_a_AAt(w, dshapex, pelmat_comp); // (grad u^{n+1},v)
+
+      for (int i = 0; i < dim; i++)
+      {
+         elmat.AddMatrix(.5,elmat_comp, i * nd, i * nd);
+         elmat.AddMatrix(-.5,elmat_comp_T, i * nd, i * nd);
+         elmat.AddMatrix(pelmat_comp, i * nd, i * nd);
+      }
+   }
+}
+
 }

fem/nonlininteg.hpp

@@ -451,6 +451,50 @@ public:
                                     DenseMatrix &elmat);
 };
 
+
+/** This class is used to assemble the first block of the Navier-Stokes stokes equation
+ * using the convective form of the nonlinear term arising in the Navier-Stokes equations \f \nu (\nabla u, \nabla v) + $(u \cdot \nabla v, w )\f$ */
+class ConvectiveNSENLFIntegrator :
+   public VectorConvectionNLFIntegrator
+{
+private:
+   Coefficient *Q{};
+   DenseMatrix dshape, dshapex, EF, gradEF, ELV, elmat_comp, pelmat_comp;
+   Vector shape;
+
+public:
+   ConvectiveNSENLFIntegrator(Coefficient &q): Q(&q) { }
+
+   ConvectiveNSENLFIntegrator() = default;
+
+   virtual void AssembleElementGrad(const FiniteElement &el,
+                                    ElementTransformation &trans,
+                                    const Vector &elfun,
+                                    DenseMatrix &elmat);
+};
+
+
+/** This class is used to assemble the first block of the Navier-Stokes stokes equation
+ * using the skew-symmetric form of the nonlinear term arising in the Navier-Stokes equations \f \nu (\nabla u, \nabla v) + $.5*(u \cdot \nabla v, w ) - .5*(u \cdot \nabla w, v )\f$ */
+class SkewSymmetricNSENLFIntegrator :
+   public VectorConvectionNLFIntegrator
+{
+private:
+   Coefficient *Q{};
+   DenseMatrix dshape, dshapex, EF, gradEF, ELV, elmat_comp, pelmat_comp;
+   Vector shape;
+
+public:
+   SkewSymmetricNSENLFIntegrator(Coefficient &q): Q(&q) { }
+
+   SkewSymmetricNSENLFIntegrator() = default;
+
+   virtual void AssembleElementGrad(const FiniteElement &el,
+                                    ElementTransformation &trans,
+                                    const Vector &elfun,
+                                    DenseMatrix &elmat);
+};
+
 }
 
 #endif