Various framework modifications related to finite deformation plasitcity. Added project method in SIMbase and NonLinSIM.

git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@871 e10b68d5-8a6e-419e-a041-bce267b0401d

src/ASM/ASMs2D.C

@@ -943,7 +943,7 @@ bool ASMs2D::integrate (Integrand& integrand,
 	}
 
       // Finalize the element quantities
-      if (!integrand.finalizeElement(elmInt))
+      if (!integrand.finalizeElement(elmInt,time))
 	return false;
 
       // Assembly of global system integral

src/ASM/ASMs2DLag.C

@@ -241,7 +241,7 @@ bool ASMs2DLag::integrate (Integrand& integrand,
 	}
 
       // Finalize the element quantities
-      if (!integrand.finalizeElement(elmInt))
+      if (!integrand.finalizeElement(elmInt,time))
 	return false;
 
       // Assembly of global system integral

src/ASM/ASMs3D.C

@@ -1265,7 +1265,7 @@ bool ASMs3D::integrate (Integrand& integrand,
 	    }
 
 	// Finalize the element quantities
-	if (!integrand.finalizeElement(elmInt))
+	if (!integrand.finalizeElement(elmInt,time))
 	  return false;
 
 	// Assembly of global system integral

src/ASM/ASMs3DLag.C

@@ -259,7 +259,7 @@ bool ASMs3DLag::integrate (Integrand& integrand,
 	    }
 
 	// Finalize the element quantities
- 	if (!integrand.finalizeElement(elmInt))
+ 	if (!integrand.finalizeElement(elmInt,time))
  	  return false;
 
 	// Assembly of global system integral

src/ASM/IntegrandBase.h

@@ -59,6 +59,8 @@ public:
   //! \details This method is invoked once before starting the numerical
   //! integration over the entire spatial domain.
   virtual void initIntegration(const TimeDomain&) {}
+  //! \brief Initializes the integrand for a new result point loop.
+  virtual void initResultPoints() {}
 
   //! \brief Initializes current element for numerical integration.
   //! \param[in] MNPC Matrix of nodal point correspondance for current element
@@ -125,7 +127,10 @@ public:
   //! It can also be used to implement multiple integration point loops within
   //! the same element, provided the necessary integration point values are
   //! stored internally in the object during the first integration loop.
-  virtual bool finalizeElement(LocalIntegral*&) { return true; }
+  virtual bool finalizeElement(LocalIntegral*&, const TimeDomain&)
+  {
+    return true;
+  }
 
   //! \brief Evaluates the integrand at an interior point.
   //! \param elmInt The local integral object to receive the contributions
@@ -350,7 +355,7 @@ public:
   //! \param[out] s The solution field value at current point
   //! \param[in] asol The analytical solution field (scalar field)
   //! \param[in] X Cartesian coordinates of current point
-  virtual bool evalPrimSol(real& s, const RealFunc& asol, const Vec3& X) const
+  virtual bool evalPrimSol(double& s, const RealFunc& asol, const Vec3& X) const
   {
     std::cerr <<" *** Integrand::evalPrimSol not implemented"<< std::endl;
     return false;

src/ASM/TimeDomain.h

@@ -27,7 +27,7 @@ struct TimeDomain
   bool first; //!< If \e true, this is the first load/time step
 
   //! \brief Default constructor.
-  TimeDomain() { t = dt = 0.0; it = 0; first = true; }
+  TimeDomain(int i = 0, bool f = true) : it(i), first(f) { t = dt = 0.0; }
 };
 
 #endif

src/Integrands/Elasticity.C

@@ -143,6 +143,9 @@ void Elasticity::setTraction (TractionFunc* tf)
 
 bool Elasticity::initElement (const std::vector<int>& MNPC)
 {
+  if (myMats)
+    myMats->withLHS = true;
+
   const size_t nen = MNPC.size();
 
   if (eKm) eKm->resize(nsd*nen,nsd*nen,true);
@@ -153,39 +156,39 @@ bool Elasticity::initElement (const std::vector<int>& MNPC)
   // Extract the current primary solution and store in the array *eV
   int ierr = 0;
   if (eV && !primsol.front().empty())
-    if ((ierr = utl::gather(MNPC,nsd,primsol.front(),*eV)))
-      std::cerr <<" *** Elasticity::initElement: Detected "
-		<< ierr <<" node numbers out of range."<< std::endl;
+    ierr = utl::gather(MNPC,nsd,primsol.front(),*eV);
 
   // If previous solutions are required, they are stored in the
   // (primsol.size()-1) last entries in myMats->b
   int j = 1+myMats->b.size()-primsol.size();
   for (size_t i = 1; i < primsol.size() && ierr == 0; i++, j++)
-    ierr = utl::gather(MNPC,nsd,primsol[i],myMats->b[j]);
+    if (!primsol[i].empty())
+      ierr = utl::gather(MNPC,nsd,primsol[i],myMats->b[j]);
 
-  if (myMats)
-    myMats->withLHS = true;
+  if (ierr == 0) return true;
 
-  return ierr == 0;
+  std::cerr <<" *** Elasticity::initElement: Detected "
+	    << ierr <<" node numbers out of range."<< std::endl;
+  return false;
 }
 
 
 bool Elasticity::initElementBou (const std::vector<int>& MNPC)
 {
-  const size_t nen = MNPC.size();
-
-  if (eS) eS->resize(nsd*nen,true);
-
-  int ierr = 0;
-  if (eV && !primsol.front().empty())
-    if ((ierr = utl::gather(MNPC,nsd,primsol.front(),*eV)))
-      std::cerr <<" *** Elasticity::initElement: Detected "
-                << ierr <<" node numbers out of range."<< std::endl;
-
   if (myMats)
     myMats->withLHS = false;
 
-  return ierr == 0;
+  if (eS) eS->resize(nsd*MNPC.size(),true);
+
+  int ierr = 0;
+  if (eV && !primsol.front().empty())
+    ierr = utl::gather(MNPC,nsd,primsol.front(),*eV);
+
+  if (ierr == 0) return true;
+
+  std::cerr <<" *** Elasticity::initElement: Detected "
+	    << ierr <<" node numbers out of range."<< std::endl;
+  return false;
 }
 
 

src/Integrands/MaterialBase.h

@@ -39,6 +39,11 @@ public:
   //! \brief Prints out material parameters to the given output stream.
   virtual void print(std::ostream&) const {}
 
+  //! \brief Initializes the material model for a new integration loop.
+  virtual void initIntegration(const TimeDomain&) {}
+  //! \brief Initializes the material model for a new result point loop.
+  virtual void initResultPoints() {}
+
   //! \brief Evaluates the mass density at current point.
   virtual double getMassDensity(const Vec3&) const { return 0.0; }
 

src/SIM/NonLinSIM.C

@@ -159,9 +159,12 @@ bool NonLinSIM::solveStep (SolvePrm& param, SIM::SolutionMode mode)
 	if (!this->updateConfiguration(param))
 	  return false;
 
-	param.time.first = false;
 	model->setMode(SIM::RECOVERY);
-	return this->solutionNorms(param.time);
+	if (!this->solutionNorms(param.time))
+	  return false;
+
+	param.time.first = false;
+	return true;
 
       case DIVERGED:
 	return false;
@@ -421,3 +424,9 @@ void NonLinSIM::dumpStep (int iStep, double time, std::ostream& os,
 
   model->dumpPrimSol(solution.front(),os,withID);
 }
+
+
+bool NonLinSIM::project ()
+{
+  return model->project(solution.front());
+}

src/SIM/NonLinSIM.h

@@ -107,6 +107,13 @@ public:
   //! \brief Returns a const reference to current solution vector.
   const Vector& getSolution() const { return solution.front(); }
 
+  //! \brief Projects the secondary solution onto the spline basis.
+  //! \details The secondary solution, defined through the Integrand class,
+  //! is projected onto the spline basis to obtain the control point values
+  //! of the secondary solution. These values then overwrite the current
+  //! primary solution vector.
+  bool project();
+
 protected:
   //! \brief Convergence status enum.
   enum ConvStatus { NONE, CONVERGED, DIVERGED };

src/SIM/SIMbase.C

@@ -610,6 +610,8 @@ bool SIMbase::solutionNorms (const TimeDomain& time, const Vectors& psol,
 
   PROFILE1("Norm integration");
 
+  myProblem->initIntegration(time);
+
   size_t nNorms = norm->getNoFields() + norm->getNoSecFields();
   const Vector& primsol = psol.front();
 
@@ -919,6 +921,9 @@ bool SIMbase::writeGlvS (const Vector& psol,
   else if (!myVtf)
     return false;
 
+  if (!psolOnly)
+    myProblem->initResultPoints();
+
   Matrix field;
   size_t i, j;
   int geomID = 0, idBlock = 10;
@@ -1283,3 +1288,18 @@ bool SIMbase::finalizeAssembly (bool newLHSmatrix)
 
   return true;
 }
+
+
+bool SIMbase::project (Vector& psol)
+{
+  Matrix values;
+  Vector psol2(psol.size());
+  for (size_t i = 0; i < myModel.size(); ++i) {
+    myModel[i]->extractNodeVec(psol,myProblem->getSolution());
+    myModel[i]->evalSolution(values,*myProblem,NULL);
+    myModel[i]->injectNodeVec(values,psol2);
+  }
+  psol = psol2;
+
+  return true;
+}

src/SIM/SIMbase.h

@@ -213,6 +213,13 @@ public:
 		   int nev, int ncv, int iop, double shift,
 		   size_t iA = 0, size_t iB = 1);
 
+  //! \brief Project the secondary solution associated with \a psol.
+  //! \param psol Control point values of primary(in)/secondary(out) solution
+  //!
+  //! \details The secondary solution, defined through the Integrand class,
+  //! corresponding to the primary solution \a psol is projected onto the
+  //! spline basis to obtain the control point values of the secondary solution.
+  bool project(Vector& psol);
 
   // Post-processing methods
   // =======================

src/Utility/Tensor.C

@@ -346,7 +346,7 @@ SymmTensor::SymmTensor (const SymmTensor& T) : Tensor(0)
 
 SymmTensor& SymmTensor::transform (const Tensor& T)
 {
-  if (T.symmetric() || T.dim() != n) return *this;
+  if (T.symmetric() || T.dim() < n) return *this;
 
   real S11, S12, S13, S21, S22, S23, S31, S32, S33;
   switch (n) {