resolve conflict

examples/makefile

@@ -72,7 +72,8 @@ SUBDIRS_TPRINT = $(addsuffix /test-print,$(SUBDIRS))
 
 # Replace the default implicit rule for *.cpp files
 %: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
-	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS) -I/Users/mittal3/local/include/.
+
 
 all: $(EXAMPLES) $(SUBDIRS_ALL)
 

examples/test_integral.cpp

@@ -0,0 +1,241 @@
+//                                MFEM Example 15
+//
+
+#include "../mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "../../algoim/src/algoim_quad.hpp"
+#include "blitz/tinyvec2.h"
+#include "../linalg/ttensor.hpp"
+
+using namespace std;
+using namespace mfem;
+using namespace Algoim;
+
+template<typename T>
+double GetTValue(T x)
+{
+   return (double)x;
+};
+
+template<int N>
+double GetTValue(Algoim::Interval<N> x)
+{
+   return x.alpha;
+};
+
+template<typename T>
+void SetTValue(T &x, double xv)
+{
+   x = (T)xv;
+};
+
+template<int N>
+void SetTValue(Algoim::Interval<N> &x, double xv)
+{
+   x.alpha = xv;
+};
+
+#define radius 0.6
+#define lstype 3
+template<int N, typename T>
+T ReturnFuncValue(blitz::TinyVector<T,N> &xc)
+{
+   T fx;
+   if (lstype == 1)
+   {
+      fx = 1*(xc(0)*xc(0) + xc(1)*xc(1) - radius*radius);
+   }
+   else if (lstype == 2)
+   {
+      double a1 = 20., a2 = 2., a3 = 3.;
+      T yv = a1*(xc(1)-0.5),
+        xv = a2*sin(a3*(xc(0)-0.5)*M_PI);
+      fx = tanh(yv + xv + 1);
+   }
+   else if (lstype == 3)
+   {
+      const int num_circ = 3;
+      double rad[num_circ] = {0.3, 0.15, 0.2};
+      double c[num_circ][2] = { {0.6, 0.6}, {0.3, 0.3}, {0.25, 0.75} };
+
+      const double xv = xc(0), yv = xc(1);
+
+      // circle 0
+      double r0 = (xv-c[0][0])*(xv-c[0][0]) + (yv-c[0][1])*(yv-c[0][1]);
+      r0 = (r0 > 0) ? std::sqrt(r0) : 0.0;
+      if (r0 <= 0.2) { return -1.0; }
+
+      for (int i = 0; i < num_circ; i++)
+      {
+         double r = (xv-c[i][0])*(xv-c[i][0]) + (yv-c[i][1])*(yv-c[i][1]);
+         r = (r > 0) ? std::sqrt(r) : 0.0;
+         if (r <= rad[i]) { return 1.0; }
+      }
+
+      // rectangle 1
+      if (0.7 <= xv && xv <= 0.8 && 0.1 <= yv && yv <= 0.8) { return 1.0; }
+
+      // rectangle 2
+      if (0.3 <= xv && xv <= 0.8 && 0.15 <= yv && yv <= 0.2) { return 1.0; }
+      return -1.0;
+   }
+   else
+   {
+      fx = (T)(0.);
+   }
+   return fx;
+};
+
+template<int N, typename T>
+blitz::TinyVector<T,N> ReturnFuncGradient(blitz::TinyVector<T,N> &xc,
+                                          DenseMatrix &J)
+{
+   blitz::TinyVector<T,N> dfx;
+   if (lstype == 1)
+   {
+      dfx = blitz::TinyVector<T,N>(2.0*J(0, 0)*xc(0) + 2.0*J(1, 0)*xc(1),
+                                   2.0*J(0, 1)*xc(0) + 2.0*J(1, 1)*xc(1));
+   }
+   else if (lstype == 2)
+   {
+      double a1 = 20., a2 = 2., a3 = 3.;
+      T yv = a1*(xc(1)-0.5),
+        xv = a2*sin(a3*(xc(0)-0.5)*M_PI),
+        scale = sech(yv+xv+1);
+      dfx  = blitz::TinyVector<T,N>(scale*xv*a3*J(0, 0) + scale*a1*J(1, 0),
+                                    scale*xv*a3*J(0, 1) + scale*a1*J(1, 1));
+   }
+   else
+   {
+      dfx = blitz::TinyVector<T,N>((T)(0.), (T)(0.));
+   }
+   return dfx;
+};
+
+double levelset(const Vector &x)
+{
+   blitz::TinyVector<double, 2> xc;
+   for (int i = 0; i < x.Size(); i++) { xc(i) = x(i); }
+   double val = ReturnFuncValue(xc);
+   return val;
+}
+
+
+template<int N>
+struct AnalyticalLevelSet
+{
+private:
+   ElementTransformation *Tr;
+   GridFunction *Gf;
+public:
+   AnalyticalLevelSet(ElementTransformation &Tr_,
+                      GridFunction &Gf_) : Tr(&Tr_), Gf(&Gf_) { }
+
+   template<typename T>
+   T operator() (blitz::TinyVector<T,N>& x) const
+   {
+      return Gf->GetTValue(Tr->ElementNo, x);
+   }
+
+   template<typename T>
+   blitz::TinyVector<T,N> grad(blitz::TinyVector<T,N>& x) const
+   {
+      blitz::TinyVector<T,N> dfx;
+      Gf->GetTGradient(Tr->ElementNo, x, dfx);
+      return dfx;
+   }
+};
+
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/inline-quad.mesh";
+   int order = 2;
+   int rs_levels         = 0;
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&rs_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   Mesh mesh_(mesh_file, 1, 1, false);
+   for (int lev = 0; lev < rs_levels; lev++) { mesh_.UniformRefinement(); }
+
+   H1_FECollection h1fec_(order, mesh_.Dimension());
+   FiniteElementSpace h1fes_(&mesh_, &h1fec_);
+   GridFunction x0(&h1fes_);
+
+   FunctionCoefficient ind(levelset);
+   x0.ProjectCoefficient(ind);
+   if (visualization)
+   {
+      osockstream sock(19916, "localhost");
+      sock << "solution\n";
+      mesh_.Print(sock);
+      x0.Save(sock);
+      sock.send();
+      sock << "window_title 'Level set'\n"
+           << "window_geometry "
+           << 1200 << " " << 0 << " " << 600 << " " << 600 << "\n"
+           << "keys jRmclA" << endl;
+   }
+
+   double area = 0.0;
+   ElementTransformation *Tr = NULL;
+   IntegrationPoint *ip = new IntegrationPoint();
+   Vector el_area(mesh_.GetNE());
+   el_area = 0.0;
+   ofstream myfile;
+   myfile.open ("qpts.out");
+   for (int e = 0; e < mesh_.GetNE(); e++)
+   {
+      Tr = mesh_.GetElementTransformation(e);
+      AnalyticalLevelSet<2> phi(*Tr, x0);
+      auto q = Algoim::quadGen<2>(phi, Algoim::BoundingBox<double,2>(0.0, 1.0), -2,
+                                  -1, order);
+      double elsum = 0.0;
+      for (const auto& pt : q.nodes)
+      {
+         ip->Set2(pt.x(0), pt.x(1));
+         Tr->SetIntPoint(ip);
+         Vector xtm(2);
+         Tr->Transform(*ip, xtm);
+         area += Tr->Weight() * pt.w;
+         elsum += Tr->Weight() * pt.w;
+         myfile << xtm(0) << " " << xtm(1) << endl; //write quadrature points to file.
+      }
+      el_area(e) = elsum;
+   }
+   myfile.close();
+   double exact_area;
+   if (lstype == 1)
+   {
+      exact_area = M_PI*radius*radius/4;
+   }
+   else if (lstype == 2)
+   {
+      exact_area = 0.45;
+   }
+   std::cout << " Location of integration points output in qpts.out.\n";
+   cout << " Numerical area: " <<  std::setprecision(5) << area << endl;
+   cout << " Exact area:     " <<  std::setprecision(5) << exact_area  << endl;
+   cout << " Error:          " <<  std::setprecision(5) << std::fabs(
+           area-exact_area) << endl;
+   //   el_area.Print();
+
+
+   return 0;
+}

fem/fe.hpp

@@ -17,6 +17,8 @@
 #include "../linalg/linalg.hpp"
 #include "intrules.hpp"
 #include "geom.hpp"
+#include "blitz/tinyvec2.h"
+#include "./../../algoim/src/algoim_interval.hpp"
 
 #include <map>
 
@@ -1987,6 +1989,71 @@ public:
       void Eval(const double x, Vector &u) const;
       void Eval(const double x, Vector &u, Vector &d) const;
       void Eval(const double x, Vector &u, Vector &d, Vector &d2) const;
+      template<int N>
+      void EvalT(const Algoim::Interval<N> &y,
+                 std::vector<Algoim::Interval<N>> &u) const
+      {
+         int i, k, p = x.Size() - 1;
+         for (i = 0; i <= p; i++)
+         {
+            u[i] = 0.0;
+         }
+         Algoim::Interval<N> lk(1.0);
+
+         for (i = 0; i <= p; i++)
+         {
+            lk = 1.0;
+            for (k = 0; k <= p; k++)
+            {
+               if (k == i) { continue; }
+               lk *= y - x(k);
+            }
+            u[i] = lk * w(i);
+         }
+      }
+
+      template<int N>
+      void EvalT(const Algoim::Interval<N> &y,
+                 std::vector<Algoim::Interval<N>> &u,
+                 std::vector<Algoim::Interval<N>> &d) const
+      {
+         int i, k, m, p = x.Size() - 1;
+         for (i = 0; i <= p; i++)
+         {
+            u[i] = 0.0;
+            d[i] = 0.0;
+         }
+         Algoim::Interval<N> lk(1.0);
+
+         for (i = 0; i <= p; i++)
+         {
+            lk = 1.0;
+            for (k = 0; k <= p; k++)
+            {
+               if (k == i) { continue; }
+               lk *= y - x(k);
+            }
+            u[i] = lk * w(i);
+         }
+
+         for (i = 0; i <= p; i++)
+         {
+            lk = 0.0;
+            for (k = 0; k <= p; k++)
+            {
+               if (k == i) { continue; }
+               Algoim::Interval<N> prod(1.0);
+               for (m = 0; m <= p; m++)
+               {
+                  if (k == m || m == i) { continue; }
+                  prod *= y - x(m);
+               }
+               lk += prod;
+            }
+            d[i] = lk * w(i);
+         }
+      }
+
       /// Evaluate the "integrated" basis, which is given by the negative
       /// partial sum of the corresponding closed basis derivatives. The closed
       /// basis derivatives are given by @a d, and the result is stored in @a i.
@@ -2290,6 +2357,50 @@ public:
    virtual void CalcHessian(const IntegrationPoint &ip,
                             DenseMatrix &Hessian) const;
    virtual void ProjectDelta(int vertex, Vector &dofs) const;
+   template<int N>
+   void CalcTShape(const blitz::TinyVector<Algoim::Interval<N>,N> &ip,
+                   std::vector<Algoim::Interval<N>> &shape) const
+   {
+      {
+         const int p = order;
+         std::vector<Algoim::Interval<N>> shape_xt(p+1), shape_yt(p+1);
+
+         basis1d.EvalT(ip(0), shape_xt);
+         basis1d.EvalT(ip(1), shape_yt);
+
+         for (int o = 0, j = 0; j <= p; j++)
+            for (int i = 0; i <= p; i++)
+            {
+               shape[dof_map[o++]] = shape_xt[i]*shape_yt[j];
+            }
+      }
+   }
+
+   template<int N>
+   void CalcTDShape(const blitz::TinyVector<Algoim::Interval<N>,N> &ip,
+                    std::vector<Algoim::Interval<N>> &dshape_x,
+                    std::vector<Algoim::Interval<N>> &dshape_y) const
+   {
+      {
+         const int p = order;
+         std::vector<Algoim::Interval<N>> shape_xt(p+1), shape_yt(p+1),
+             dshape_xt(p+1), dshape_yt(p+1);;
+
+         basis1d.EvalT(ip(0), shape_xt, dshape_xt);
+         basis1d.EvalT(ip(1), shape_yt, dshape_yt);
+
+         for (int o = 0, j = 0; j <= p; j++)
+         {
+            for (int i = 0; i <= p; i++)
+            {
+               dshape_x[dof_map[o]] = dshape_xt[i]* shape_yt[j];
+               dshape_y[dof_map[o]] =  shape_xt[i]*dshape_yt[j];
+               o++;
+            }
+         }
+      }
+   }
+
 };
 
 

fem/gridfunc.cpp

@@ -424,6 +424,18 @@ const
    return (DofVal * LocVec);
 }
 
+//template<int N>
+//Algoim::Interval<N> GridFunction::GetTValue(int i,
+//                                            blitz::TinyVector<Algoim::Interval<N>,N> &ip)
+//{
+//   return (Algoim::Interval<N>)(1.0);
+//}
+
+//double GridFunction::GetTValue(int i, blitz::TinyVector<double, 2> &ip)
+//{
+//   return (double)(1.0);
+//}
+
 void GridFunction::GetVectorValue(int i, const IntegrationPoint &ip,
                                   Vector &val) const
 {
@@ -1660,6 +1672,27 @@ void GridFunction::GetGradients(ElementTransformation &tr,
    }
 }
 
+void GridFunction::GetGradient(ElementTransformation &tr,
+                               const IntegrationPoint &ip,
+                               Vector &grad) const
+{
+   int elNo = tr.ElementNo;
+   const FiniteElement *fe = fes->GetFE(elNo);
+   MFEM_ASSERT(fe->GetMapType() == FiniteElement::VALUE, "invalid FE map type");
+   DenseMatrix dshape(fe->GetDof(), fe->GetDim());
+   Vector lval, gh(fe->GetDim()), gcol;
+   Array<int> dofs;
+   fes->GetElementDofs(elNo, dofs);
+   GetSubVector(dofs, lval);
+   grad.SetSize(fe->GetDim());
+
+   fe->CalcDShape(ip, dshape);
+   dshape.MultTranspose(lval, gh);
+   tr.SetIntPoint(&ip);
+   const DenseMatrix &Jinv = tr.InverseJacobian();
+   Jinv.MultTranspose(gh, grad);
+}
+
 void GridFunction::GetVectorGradient(
    ElementTransformation &T, DenseMatrix &grad) const
 {

fem/gridfunc.hpp

@@ -16,6 +16,8 @@
 #include "fespace.hpp"
 #include "coefficient.hpp"
 #include "bilininteg.hpp"
+#include "blitz/tinyvec2.h"
+#include "../../algoim/src/algoim_quad.hpp"
 #ifdef MFEM_USE_ADIOS2
 #include "../general/adios2stream.hpp"
 #endif
@@ -164,6 +166,98 @@ public:
    virtual double GetValue(int i, const IntegrationPoint &ip,
                            int vdim = 1) const;
 
+   template<int N, typename T>
+   T GetTValue(const int i, const blitz::TinyVector<T,N> &ipt)
+   {
+      MFEM_VERIFY(N <= 3 && N >1, "TinyVector size must be 1, 2, or 3.");
+      IntegrationPoint *ip = new IntegrationPoint();
+      if (N == 2)
+      {
+         ip->Set2(ipt(0), ipt(1));
+      }
+      else if (N == 3)
+      {
+         ip->Set3(ipt(0), ipt(1), ipt(2));
+      }
+      double value = GetValue(i, *ip);
+      return (T)value;
+   }
+
+   template<int N>
+   Algoim::Interval<N> GetTValue(const int i,
+                                 const blitz::TinyVector<Algoim::Interval<N>,N> &ip)
+   {
+      Array<int> dofs;
+      fes->GetElementDofs(i, dofs);
+      fes->DofsToVDofs(0, dofs);
+      Vector LocVec;
+      std::vector<Algoim::Interval<N>> DofValT(dofs.Size());
+      const H1_QuadrilateralElement *fe =
+         dynamic_cast<const H1_QuadrilateralElement *>(fes->GetFE(i));
+      if (fe && fe->GetMapType() == FiniteElement::VALUE)
+      {
+         fe->CalcTShape(ip, DofValT);
+      }
+      else
+      {
+         MFEM_ABORT(" Only Value type Map supported right now.");
+      }
+      GetSubVector(dofs, LocVec);
+      Algoim::Interval<N> dotp(0.0);
+      for (int i = 0; i < LocVec.Size(); i++)
+      {
+         dotp += LocVec(i)*DofValT[i];
+      }
+
+      return dotp;
+   }
+
+   template<int N, typename T>
+   void GetTGradient(const int i,
+                     const blitz::TinyVector<T,N> &ipt,
+                     blitz::TinyVector<T,N> dfx)
+   {
+      MFEM_VERIFY(N <= 3 && N >1, "TinyVector size must be 1, 2, or 3.");
+      IntegrationPoint *ip = new IntegrationPoint();
+      if (N == 2)
+      {
+         ip->Set2(ipt(0), ipt(1));
+      }
+      else if (N == 3)
+      {
+         ip->Set3(ipt(0), ipt(1), ipt(2));
+      }
+      Vector grad(N);
+      GetGradient(*fes->GetElementTransformation(i), *ip, grad);
+      for (int i = 0; i < N; i++) { dfx(i) = (T)grad(i); }
+   }
+
+   template<int N>
+   void GetTGradient(const int i,
+                     const blitz::TinyVector<Algoim::Interval<N>,N> &ip,
+                     blitz::TinyVector<Algoim::Interval<N>,N> dfx)
+   {
+      const H1_QuadrilateralElement *fe =
+         dynamic_cast<const H1_QuadrilateralElement *>(fes->GetFE(i));
+      MFEM_ASSERT(fe->GetMapType() == FiniteElement::VALUE, "invalid FE map type");
+      Vector lval;
+      Array<int> dofs;
+      fes->GetElementDofs(i, dofs);
+      GetSubVector(dofs, lval);
+
+      std::vector<Algoim::Interval<N>> dshape_x(dofs.Size()), dshape_y(dofs.Size());
+
+      fe->CalcTDShape(ip, dshape_x, dshape_y);
+      Algoim::Interval<N> u_x(0.0), u_y(0.0);
+      for (int i = 0; i < lval.Size(); i++)
+      {
+         u_x += lval(i)*dshape_x[i];
+         u_y += lval(i)*dshape_y[i];
+      }
+      dfx(0) = u_x;
+      dfx(1) = u_y;
+   }
+
    /** Return a vector value from within the given element. */
    virtual void GetVectorValue(int i, const IntegrationPoint &ip,
                                Vector &val) const;
@@ -317,6 +411,9 @@ public:
    void GetGradients(ElementTransformation &tr, const IntegrationRule &ir,
                      DenseMatrix &grad) const;
 
+   void GetGradient(ElementTransformation &tr, const IntegrationPoint &ip,
+                    Vector &grad) const;
+
    void GetGradients(const int elem, const IntegrationRule &ir,
                      DenseMatrix &grad) const
    { GetGradients(*fes->GetElementTransformation(elem), ir, grad); }

makefile

@@ -427,7 +427,7 @@ MFEM_BUILD_FLAGS = $(MFEM_PICFLAG) $(MFEM_CPPFLAGS) $(MFEM_CXXFLAGS)\
 
 # Rules for compiling all source files.
 $(OBJECT_FILES): $(BLD)%.o: $(SRC)%.cpp $(CONFIG_MK)
-	$(MFEM_CXX) $(MFEM_BUILD_FLAGS) -c $(<) -o $(@)
+	$(MFEM_CXX) $(MFEM_BUILD_FLAGS) -c $(<) -o $(@) -I/Users/mittal3/local/include/.
 
 all: examples miniapps $(TEST_DIRS)