Added: Element activation function

Author: kmokstad

src/ASM/ASMbase.C

@@ -67,6 +67,7 @@ ASMbase::ASMbase (unsigned char n_p, unsigned char n_s, unsigned char n_f)
   nel = nnod = 0;
   idx = 0;
   firstIp = 0;
+  myElActive = nullptr;
 }
 
 
@@ -85,6 +86,7 @@ ASMbase::ASMbase (const ASMbase& patch, unsigned char n_f)
   idx = patch.idx;
   firstIp = patch.firstIp;
   // Note: Properties are _not_ copied
+  myElActive = nullptr; // Element activation function is not copied
 }
 
 
@@ -123,6 +125,7 @@ ASMbase::ASMbase (const ASMbase& patch)
     std::cerr <<"  ** ASMbase copy constructor: The copied patch has"
               <<" multi-point constraints, these are not copied.\n";
 
+  myElActive = nullptr; // Element activation function is not copied
   nLag = 0; // Lagrange multipliers are not copied
 }
 
@@ -131,6 +134,7 @@ ASMbase::~ASMbase ()
 {
   for (MPC* mpc : mpcs)
     delete mpc;
+  delete myElActive;
 }
 
 
@@ -1747,10 +1751,16 @@ void ASMbase::getBoundaryElms (int lIndex, IntVec& elms,
 }
 
 
-bool ASMbase::isElementActive (int elmId) const
+bool ASMbase::isElementActive (int elmId, double time) const
 {
-  if (elmId < 1) return false;
-  if (!myActiveEls) return true;
+  if (elmId < 1)
+    return false; // element with zero extension
+
+  if (myElActive && time < (*myElActive)(elmId))
+    return false; // element not activated yet
+
+  if (!myActiveEls)
+    return true; // all elements are active
 
   return std::find(myActiveEls->begin(),
                    myActiveEls->end(),elmId) != myActiveEls->end();

src/ASM/ASMbase.h

@@ -21,11 +21,12 @@
 #include <array>
 #include <string>
 
-typedef std::vector<int>       IntVec;  //!< General integer vector
-typedef std::vector<IntVec>    IntMat;  //!< General 2D integer matrix
-typedef std::map<MPC*,int,MPCLess> MPCMap; //!< MPC-to-function code mapping
-typedef std::set<MPC*,MPCLess> MPCSet;  //!< Sorted set of MPC-equations
-typedef MPCSet::const_iterator MPCIter; //!< Iterator over an MPC-equation set
+using IntVec = std::vector<int>;    //!< General integer vector
+using IntMat = std::vector<IntVec>; //!< General 2D integer matrix
+
+using MPCMap  = std::map<MPC*,int,MPCLess>; //!< MPC-to-function code mapping
+using MPCSet  = std::set<MPC*,MPCLess>; //!< Sorted set of MPC-equations
+using MPCIter = MPCSet::const_iterator; //!< Iterator over an MPC-equation set
 
 struct TimeDomain;
 class ElementBlock;
@@ -43,9 +44,12 @@ class RealFunc;
 class VecFunc;
 class Vec3;
 class Tensor;
+
 namespace ASM { class InterfaceChecker; }
+namespace utl { template<class Arg, class Result> class Function; }
 
-typedef std::vector<ASMbase*> ASMVec; //!< Spline patch container
+using ASMVec  = std::vector<ASMbase*>;     //!< Spline patch container type
+using IntFunc = utl::Function<int,double>; //!< Real-valued integer function
 
 
 /*!
@@ -78,7 +82,7 @@ class ASMbase
     explicit BC(int n) : node(n), CX(1), CY(1), CZ(1), RX(1), RY(1), RZ(1) {}
   };
 
-  typedef std::vector<BC> BCVec; //!< Nodal boundary condition container
+  using BCVec = std::vector<BC>; //!< Nodal boundary condition container
 
 protected:
   //! \brief The constructor sets the number of space dimensions and fields.
@@ -167,6 +171,9 @@ class ASMbase
   //! It must be invoked only before SIMbase::preprocess() is invoked.
   void setNoFields(unsigned char n) { nf = n; }
 
+  //! \brief Sets the element activation function.
+  void setElementActivator(IntFunc* efunc) { myElActive = efunc; }
+
   //! \brief Sets the minimum element size for adaptive refinement.
   virtual void setMinimumSize(double) {}
   //! \brief Defines the minimum element size for adaptive refinement.
@@ -318,7 +325,7 @@ class ASMbase
   //! \brief Sets the list of active elements during assembly.
   void setActiveElements(IntVec* active) { myActiveEls = active; }
   //! \brief Returns \e true if element with global id \a elmId is active.
-  bool isElementActive(int elmId) const;
+  bool isElementActive(int elmId, double time = -1.0) const;
 
   //! \brief Returns the nodal point correspondance array for an element.
   //! \param[in] iel 1-based element index local to current patch
@@ -1024,11 +1031,12 @@ class ASMbase
   std::vector<char> myLMTypes; //!< Type of %Lagrange multiplier ('L' or 'G')
   std::set<size_t>  myLMs;     //!< Nodal indices of the %Lagrange multipliers
 
+  IntFunc* myElActive; //!< Function returning activatiation time of element
   IntVec* myActiveEls; //!< List of active elements during element assembly
   static IntVec Empty; //!< Empty integer vector used when a reference is needed
 
 protected:
-  typedef std::array<double,3> XYZ; //!< Convenience type definition
+  using XYZ = std::array<double,3>; //!< Convenience type definition
   std::map<size_t,XYZ>   myRmaster; //!< Rigid master nodal points
 };
 

src/ASM/ASMs2D.C

@@ -1755,7 +1755,8 @@ bool ASMs2D::integrate (Integrand& integrand,
       {
         int iel = groups[g][t][e];
         fe.iel = MLGE[iel];
-        if (!this->isElementActive(fe.iel)) continue; // zero-area element
+        if (!this->isElementActive(fe.iel,time.t))
+          continue; // zero-area or inactive element
 
 #ifdef SP_DEBUG
         if (dbgElm < 0 && 1+iel != -dbgElm)
@@ -2039,10 +2040,12 @@ bool ASMs2D::integrate (Integrand& integrand,
       for (size_t e = 0; e < groups[g][t].size() && ok; e++)
       {
         int iel = groups[g][t][e];
-        if (itgPts[iel].empty()) continue; // no points in this element
+        if (itgPts[iel].empty())
+          continue; // no points in this element
 
         fe.iel = MLGE[iel];
-        if (!this->isElementActive(fe.iel)) continue; // zero-area element
+        if (!this->isElementActive(fe.iel,time.t))
+          continue; // zero-area or inactive element
 
 #ifdef SP_DEBUG
         if (dbgElm < 0 && 1+iel != -dbgElm)
@@ -2218,7 +2221,8 @@ bool ASMs2D::integrate (Integrand& integrand,
       if (!hasInterfaceElms) jel = iel;
 
       fe.iel = abs(MLGE[jel]);
-      if (!this->isElementActive(fe.iel)) continue; // zero-area element
+      if (!this->isElementActive(fe.iel,time.t))
+        continue; // zero-area element
 
       short int status = iChk.hasContribution(iel,i1,i2);
       if (!status) continue; // no interface contributions for this element
@@ -2421,7 +2425,8 @@ bool ASMs2D::integrate (Integrand& integrand, int lIndex,
     for (int i1 = p1; i1 <= n1; i1++, iel++)
     {
       fe.iel = abs(MLGE[doXelms+iel-1]);
-      if (!this->isElementActive(fe.iel)) continue; // zero-area element
+      if (!this->isElementActive(fe.iel,time.t))
+        continue; // zero-area element
 
       if (!myElms.empty() && !glInt.threadSafe() &&
           std::find(myElms.begin(), myElms.end(), iel-1) == myElms.end())

src/ASM/ASMs3D.C

@@ -2110,7 +2110,8 @@ bool ASMs3D::integrate (Integrand& integrand,
       {
         int iel = groups[g][t][l];
         fe.iel = MLGE[iel];
-        if (!this->isElementActive(fe.iel)) continue; // zero-volume element
+        if (!this->isElementActive(fe.iel,time.t))
+          continue; // zero-volume or inactive element
 
 #ifdef SP_DEBUG
         if (dbgElm < 0 && 1+iel != -dbgElm)
@@ -2383,10 +2384,12 @@ bool ASMs3D::integrate (Integrand& integrand,
       for (size_t e = 0; e < groups[g][t].size() && ok; e++)
       {
         int iel = groups[g][t][e];
-        if (itgPts[iel].empty()) continue; // no points in this element
+        if (itgPts[iel].empty())
+          continue; // no integration points in this element
 
         fe.iel = MLGE[iel];
-        if (!this->isElementActive(fe.iel)) continue; // zero-volume element
+        if (!this->isElementActive(fe.iel,time.t))
+          continue; // zero-volume or inactive element
 
 #ifdef SP_DEBUG
         if (dbgElm < 0 && 1+iel != -dbgElm)
@@ -2642,7 +2645,8 @@ bool ASMs3D::integrate (Integrand& integrand, int lIndex,
       {
         int iel = threadGrp[g][t][l];
         fe.iel = abs(MLGE[doXelms+iel]);
-        if (!this->isElementActive(fe.iel)) continue; // zero-volume element
+        if (!this->isElementActive(fe.iel,time.t))
+          continue; // zero-volume or inactive element
 
 #ifdef SP_DEBUG
         if (dbgElm < 0 && 1+iel != -dbgElm)
@@ -2875,7 +2879,8 @@ bool ASMs3D::integrateEdge (Integrand& integrand, int lEdge,
       for (int i1 = p1; i1 <= n1; i1++, iel++)
       {
         fe.iel = MLGE[iel-1];
-        if (!this->isElementActive(fe.iel)) continue; // zero-volume element
+        if (!this->isElementActive(fe.iel,time.t))
+          continue; // zero-volume or inactive element
 
         if (!myElms.empty() &&
             std::find(myElms.begin(), myElms.end(), fe.iel-1) == myElms.end())

src/SIM/SIMinput.C

@@ -349,6 +349,28 @@ bool SIMinput::parseGeometryTag (const tinyxml2::XMLElement* elem)
     }
   }
 
+  else if (!strcasecmp(elem->Value(),"activation") && elem->FirstChild())
+  {
+    int patch = 0;
+    std::string type("expression");
+    utl::getAttribute(elem,"patch",patch);
+    utl::getAttribute(elem,"type",type,true);
+    if (patch < 1 || patch > nGlPatches)
+    {
+      std::cerr <<" *** SIMinput::parse: Invalid patch index "
+                << patch <<"."<< std::endl;
+      return false;
+    }
+    ASMbase* pch = this->getPatch(patch,true);
+    if (pch)
+    {
+      const char* funcdef = elem->FirstChild()->Value();
+      IFEM::cout <<"\tActivation function for P"<< patch
+                 <<": "<< funcdef << std::endl;
+      pch->setElementActivator(utl::parseIntFunc(funcdef,type));
+    }
+  }
+
   return true;
 }
 
@@ -685,8 +707,11 @@ bool SIMinput::parseICTag (const tinyxml2::XMLElement* elem)
 bool SIMinput::parseLinSolTag (const tinyxml2::XMLElement* elem)
 {
   if (!strcasecmp(elem->Value(),"class"))
+  {
     if (elem->FirstChild())
       opt.setLinearSolver(elem->FirstChild()->Value());
+    utl::getAttribute(elem,"fixzeros",fixZeros);
+  }
 
   return true;
 }