Added: Calculation of updated element coordinates

when the integration type contains UPDATED_NODES

src/ASM/ASMbase.C

@@ -1176,6 +1176,54 @@ bool ASMbase::getSolution (Matrix& sField, const Vector& locSol,
 }
 
 
+/*!
+  This method calculates the updated (control point) coordinates of current
+  element by adding the current displacement vector stored as the first vector
+  in \a eVec to the nodal coordinates of the reference configuration.
+  The result is stored as the last vector in \a eVec.
+*/
+
+bool ASMbase::deformedConfig (const RealArray& Xnod,
+                              Vectors& eVec, bool force2nd) const
+{
+  if (eVec.empty() || eVec.front().empty())
+    return true; // No primary solution yet, silently OK
+
+  const Vector& eV = eVec.front();
+  size_t nen = Xnod.size()/nsd;
+  size_t npv = nen > 0 ? eV.size()/nen : 0;
+  if (npv < nsd || npv*nen != eV.size() || nsd*nen != Xnod.size())
+  {
+    std::cerr <<" *** ASMbase::deformedConfig: Inconsistent size of "
+              <<" nodal coordinate and displacement arrays "
+              << eV.size() <<", "<< Xnod.size() << std::endl;
+    return false;
+  }
+
+  Vector Xdef(Xnod.data(),Xnod.size());
+  if (npv == nsd)
+    Xdef.add(eV);
+  else
+  {
+    size_t i, j;
+    for (i = j = 0; i < Xdef.size(); i++)
+    {
+      Xdef[i] += eV[j+i%nsd];
+      if ((i+1)%nsd == 0) j += npv;
+    }
+  }
+
+  if (force2nd && eVec.size() >= 2)
+    eVec[1].swap(Xdef);
+  else
+    eVec.push_back(Xdef);
+#if SP_DEBUG > 2
+  std::cout <<"Element solution vector "<< eVec.size() << eVec.back();
+#endif
+  return true;
+}
+
+
 int ASMbase::searchCtrlPt (RealArray::const_iterator cit,
                            RealArray::const_iterator end,
                            const Vec3& X, int dimension, double tol) const

src/ASM/ASMbase.h

@@ -724,18 +724,25 @@ public:
   //! \param[in] all Returns \e true only if all DOFs are fixed
   bool isFixed(int node, int dof, bool all = false) const;
 
+protected:
+  //! \brief Calculated the deformed configuration for current element.
+  //! \param[in] Xnod Array of nodal point coordinates for current element
+  //! \param eVec Current element solution vectors
+  //! \param[in] force2nd If \e true, put updated coordinates as the 2nd vector
+  bool deformedConfig(const RealArray& Xnod, Vectors& eVec,
+                      bool force2nd = false) const;
+
+  //! \brief Collapses the given two nodes into one.
+  //! \details The global node number of the node with the highest number
+  //! is changed into the number of the other node.
+  static bool collapseNodes(ASMbase& pch1, int node1, ASMbase& pch2, int node2);
+
 private:
   //! \brief Recursive method used by \a resolveMPCchains.
   //! \param[in] allMPCs All multi-point constraint equations in the model
   //! \param mpc Pointer to the multi-point constraint equation to resolve
   static bool resolveMPCchain(const MPCSet& allMPCs, MPC* mpc);
 
-protected:
-  //! \brief Collapses the given two nodes into one.
-  //! \details The global node number of the node with the highest number
-  //! is changed into the number of the other node.
-  static bool collapseNodes(ASMbase& pch1, int node1, ASMbase& pch2, int node2);
-
 public:
   static bool fixHomogeneousDirichlet; //!< If \e true, pre-eliminate fixed DOFs
 

src/ASM/ASMs2D.C

@@ -1649,6 +1649,15 @@ bool ASMs2D::integrate (Integrand& integrand,
           break;
         }
 
+        if (integrand.getIntegrandType() & Integrand::UPDATED_NODES)
+          if (!time.first || time.it > 0)
+            if (!this->deformedConfig(Xnod,A->vec))
+            {
+              A->destruct();
+              ok = false;
+              break;
+            }
+
         if (xr)
         {
           // --- Selective reduced integration loop ----------------------------
@@ -1884,6 +1893,15 @@ bool ASMs2D::integrate (Integrand& integrand,
           break;
         }
 
+        if (integrand.getIntegrandType() & Integrand::UPDATED_NODES)
+          if (!time.first || time.it > 0)
+            if (!this->deformedConfig(Xnod,A->vec))
+            {
+              A->destruct();
+              ok = false;
+              break;
+            }
+
 
         // --- Integration loop over all quadrature points in this element -----
 
@@ -2630,6 +2648,12 @@ bool ASMs2D::evalSolution (Matrix& sField, const IntegrandBase& integrand,
   // Fetch nodal (control point) coordinates
   Matrix Xnod, Xtmp;
   this->getNodalCoordinates(Xnod);
+  if (integrand.getIntegrandType() & Integrand::UPDATED_NODES)
+  {
+    Vectors& eV = const_cast<IntegrandBase&>(integrand).getSolutions();
+    if (!this->deformedConfig(Xnod,eV,true))
+      return false;
+  }
 
   FiniteElement fe(p1*p2,firstIp);
   fe.p = p1 - 1;

src/ASM/Integrand.h

@@ -146,14 +146,15 @@ public:
     STANDARD           = 0, //!< Default integrand type (first derivatives only)
     NO_DERIVATIVES     = 1, //!< Integrand don't want any derivatives
     SECOND_DERIVATIVES = 2, //!< Integrand wants second derivatives
-    AVERAGE            = 4, //!< Integrand wants basis function averages
-    ELEMENT_CORNERS    = 8, //!< Integrand wants element corner coordinates
-    ELEMENT_CENTER    = 16, //!< Integrand wants element center coordinates
-    G_MATRIX          = 32, //!< Integrand wants the G matrix
-    NODAL_ROTATIONS   = 64, //!< Integrand wants nodal rotation tensors
-    XO_ELEMENTS      = 128, //!< Integrand is defined on extraordinary elements
-    INTERFACE_TERMS  = 256, //!< Integrand has element interface terms
-    NORMAL_DERIVS    = 512  //!< Integrand p'th order normal derivatives
+    AVERAGE         = 1<<2, //!< Integrand wants basis function averages
+    ELEMENT_CORNERS = 1<<3, //!< Integrand wants element corner coordinates
+    ELEMENT_CENTER  = 1<<4, //!< Integrand wants element center coordinates
+    G_MATRIX        = 1<<5, //!< Integrand wants the G matrix
+    NODAL_ROTATIONS = 1<<6, //!< Integrand wants nodal rotation tensors
+    XO_ELEMENTS     = 1<<7, //!< Integrand is defined on extraordinary elements
+    INTERFACE_TERMS = 1<<8, //!< Integrand has element interface terms
+    NORMAL_DERIVS   = 1<<9, //!< Integrand p'th order normal derivatives
+    UPDATED_NODES   = 1<<10 //!< Integrand wants updated nodal coordinates
   };
 
   //! \brief Defines which FE quantities are needed by the integrand.

src/ASM/IntegrandBase.h

@@ -41,7 +41,8 @@ class IntegrandBase : public Integrand
 {
 protected:
   //! \brief The default constructor is protected to allow sub-classes only.
-  explicit IntegrandBase(unsigned short int n = 0) : nsd(n), npv(1), m_mode(SIM::INIT) {}
+  explicit IntegrandBase(unsigned short int n = 0) : nsd(n), npv(1),
+                                                     m_mode(SIM::INIT) {}
 
 public:
   //! \brief Empty destructor.

src/ASM/LR/ASMu2D.C

@@ -1105,6 +1105,15 @@ bool ASMu2D::integrate (Integrand& integrand,
         continue;
       }
 
+      if (integrand.getIntegrandType() & Integrand::UPDATED_NODES)
+        if (!time.first || time.it > 0)
+          if (!this->deformedConfig(Xnod,A->vec))
+          {
+            A->destruct();
+            ok = false;
+            continue;
+          }
+
       if (xr)
       {
         // --- Selective reduced integration loop ------------------------------
@@ -1326,6 +1335,15 @@ bool ASMu2D::integrate (Integrand& integrand,
         continue;
       }
 
+      if (integrand.getIntegrandType() & Integrand::UPDATED_NODES)
+        if (!time.first || time.it > 0)
+          if (!this->deformedConfig(Xnod,A->vec))
+          {
+            A->destruct();
+            ok = false;
+            break;
+          }
+
       // --- Integration loop over all quadrature points in this element -------
 
       size_t jp = MPitg[iel-1]; // Patch-wise integration point counter