Added support for the <neumann> XML tag on SIMbase level within the

<boundaryconditions> tag. The new tag takes all the old formats for
specifying surface pressure variations in addition to expressions.
This makes the <constantpressure> and <linearpressure> tags in the
elasticity solver obsolete and have therefore been removed.

Function expressions are now also available for Dirichlet conditions.

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

Apps/FiniteDefElasticity/NonlinearElasticityTL.C

@@ -187,7 +187,7 @@ bool NonlinearElasticityTL::evalBou (LocalIntegral& elmInt,
 				     const FiniteElement& fe,
 				     const Vec3& X, const Vec3& normal) const
 {
-  if (!tracFld)
+  if (!tracFld && !fluxFld)
   {
     std::cerr <<" *** NonlinearElasticityTL::evalBou: No tractions."
 	      << std::endl;
@@ -203,15 +203,18 @@ bool NonlinearElasticityTL::evalBou (LocalIntegral& elmInt,
   ElmMats& elMat = static_cast<ElmMats&>(elmInt);
 
   // Evaluate the surface traction
-  Vec3 T = (*tracFld)(X,normal);
+  Vec3 T = this->getTraction(X,normal);
 
   // Store traction value for visualization
-  if (fe.iGP < tracVal.size())
+  if (fe.iGP < tracVal.size() && !T.isZero())
-    if (!T.isZero()) tracVal[fe.iGP] = std::make_pair(X,T);
+  {
+    tracVal[fe.iGP].first = X;
+    tracVal[fe.iGP].second += T;
+  }
 
   // Check for with-rotated pressure load
   unsigned short int i, j;
-  if (tracFld->isNormalPressure())
+  if (tracFld && tracFld->isNormalPressure())
   {
     // Compute the deformation gradient, F
     Matrix B;

Apps/FiniteDefElasticity/NonlinearElasticityUL.C

@@ -258,7 +258,7 @@ bool NonlinearElasticityUL::evalBou (LocalIntegral& elmInt,
 				     const FiniteElement& fe,
 				     const Vec3& X, const Vec3& normal) const
 {
-  if (!tracFld)
+  if (!tracFld && !fluxFld)
   {
     std::cerr <<" *** NonlinearElasticityUL::evalBou: No tractions."
 	      << std::endl;
@@ -272,11 +272,14 @@ bool NonlinearElasticityUL::evalBou (LocalIntegral& elmInt,
   }
 
   // Evaluate the surface traction
-  Vec3 T = (*tracFld)(X,normal);
+  Vec3 T = this->getTraction(X,normal);
 
   // Store traction value for visualization
-  if (fe.iGP < tracVal.size())
+  if (fe.iGP < tracVal.size() && !T.isZero())
-    if (!T.isZero()) tracVal[fe.iGP] = std::make_pair(X,T);
+  {
+    tracVal[fe.iGP].first = X;
+    tracVal[fe.iGP].second += T;
+  }
 
   // Axi-symmetric integration point volume; 2*pi*r*|J|*w
   double detJW = axiSymmetry ? 2.0*M_PI*X.x*fe.detJxW : fe.detJxW;
@@ -290,7 +293,7 @@ bool NonlinearElasticityUL::evalBou (LocalIntegral& elmInt,
       return false;
 
     // Check for with-rotated pressure load
-    if (tracFld->isNormalPressure())
+    if (tracFld && tracFld->isNormalPressure())
     {
       // Compute its inverse and determinant, J
       double J = F.inverse();

Apps/LinearElasticity/SIMLinEl2D.C

@@ -301,26 +301,6 @@ bool SIMLinEl2D::parse (const TiXmlElement* elem)
                   << E <<" "<< nu <<" "<< rho << std::endl;
     }
 
-    else if (!strcasecmp(child->Value(),"constantpressure") ||
-             !strcasecmp(child->Value(),"linearpressure")) {
-      if (child->FirstChild() && child->FirstChild()->Value()) {
-        double p = atof(child->FirstChild()->Value());
-        int code = 0, pdir = 1;
-        utl::getAttribute(child,"code",code);
-        utl::getAttribute(child,"dir",pdir);
-        setPropertyType(code,Property::NEUMANN);
-        if (!strcasecmp(child->Value(),"linearpressure")) {
-          RealFunc* pfl = new ConstTimeFunc(new LinearFunc(p));
-          myTracs[code] = new PressureField(pfl,pdir);
-        }
-        else
-          myTracs[code] = new PressureField(p,pdir);
-        if (myPid == 0)
-          std::cout <<"\tPressure code "<< code <<" direction "<< pdir
-                    <<": "<< p << std::endl;
-      }
-    }
-
     else if (!strcasecmp(child->Value(),"anasol")) {
       std::string type;
       utl::getAttribute(child,"type",type,true);
@@ -436,12 +416,18 @@ bool SIMLinEl2D::initMaterial (size_t propInd)
 
 bool SIMLinEl2D::initNeumann (size_t propInd)
 {
-  TracFuncMap::const_iterator tit = myTracs.find(propInd);
-  if (tit == myTracs.end()) return false;
-
   Elasticity* elp = dynamic_cast<Elasticity*>(myProblem);
   if (!elp) return false;
 
-  elp->setTraction(tit->second);
+  VecFuncMap::const_iterator  vit = myVectors.find(propInd);
+  TracFuncMap::const_iterator tit = myTracs.find(propInd);
+
+  if (vit != myVectors.end())
+    elp->setTraction(vit->second);
+  else if (tit != myTracs.end())
+    elp->setTraction(tit->second);
+  else
+    return false;
+
   return true;
 }

Apps/LinearElasticity/SIMLinEl3D.C

@@ -433,26 +433,6 @@ bool SIMLinEl3D::parse (const TiXmlElement* elem)
                   << E <<" "<< nu <<" "<< rho << std::endl;
     }
 
-    else if (!strcasecmp(child->Value(),"constantpressure") ||
-             !strcasecmp(child->Value(),"linearpressure")) {
-      if (child->FirstChild() && child->FirstChild()->Value()) {
-        double p = atof(child->FirstChild()->Value());
-        int code = 0, pdir = 1;
-        utl::getAttribute(child,"code",code);
-        utl::getAttribute(child,"dir",pdir);
-        setPropertyType(code,Property::NEUMANN);
-        if (!strcasecmp(child->Value(),"linearpressure")) {
-          RealFunc* pfl = new ConstTimeFunc(new LinearFunc(p));
-          myTracs[code] = new PressureField(pfl,pdir);
-        }
-        else
-          myTracs[code] = new PressureField(p,pdir);
-        if (myPid == 0)
-          std::cout <<"\tPressure code "<< code <<" direction "<< pdir
-                    <<": "<< p << std::endl;
-      }
-    }
-
     else if (!strcasecmp(child->Value(),"anasol")) {
       std::string type;
       utl::getAttribute(child,"type",type,true);
@@ -564,12 +544,18 @@ bool SIMLinEl3D::initMaterial (size_t propInd)
 
 bool SIMLinEl3D::initNeumann (size_t propInd)
 {
-  TracFuncMap::const_iterator tit = myTracs.find(propInd);
-  if (tit == myTracs.end()) return false;
-
   Elasticity* elp = dynamic_cast<Elasticity*>(myProblem);
   if (!elp) return false;
 
-  elp->setTraction(tit->second);
+  VecFuncMap::const_iterator  vit = myVectors.find(propInd);
+  TracFuncMap::const_iterator tit = myTracs.find(propInd);
+
+  if (vit != myVectors.end())
+    elp->setTraction(vit->second);
+  else if (tit != myTracs.end())
+    elp->setTraction(tit->second);
+  else
+    return false;
+
   return true;
 }

Apps/LinearElasticity/Test/PipeJoint-Lagrange.xinp

@@ -0,0 +1,26 @@
+<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
+
+<!-- Pipe joint with shear-loaded brace.
+     Static linear-elastic analysis.
+     10-patch model, cubic Lagrange elements. !-->
+
+<simulation>
+
+  <discretization type="Lagrange"/>
+
+  <geometry>
+    <patchfile>pipe_bifurcation.g2</patchfile>
+    <nodefile>pipe_bifurcation.gno</nodefile>
+  </geometry>
+
+  <boundaryconditions>
+    <propertyfile>pipe_bifurcation.prc</propertyfile>
+    <dirichlet code="123"/>
+    <neumann code="1001" direction="1">1.0e8</neumann>
+  </boundaryconditions>
+
+  <postprocessing>
+    <vtfformat>BINARY</vtfformat>
+  </postprocessing>
+
+</simulation>

src/Integrands/Elasticity.C

@@ -41,6 +41,7 @@ Elasticity::Elasticity (unsigned short int n, bool ax) : nsd(n), axiSymmetry(ax)
   material = 0;
   locSys = 0;
   tracFld = 0;
+  fluxFld = 0;
   bodyFld = 0;
   eM = eKm = eKg = 0;
   eS = iS = 0;
@@ -173,7 +174,9 @@ LocalIntegral* Elasticity::getLocalIntegral (size_t nen, size_t,
 
 Vec3 Elasticity::getTraction (const Vec3& X, const Vec3& n) const
 {
-  if (tracFld)
+  if (fluxFld)
+    return (*fluxFld)(X);
+  else if (tracFld)
     return (*tracFld)(X,n);
   else
     return Vec3();
@@ -195,6 +198,7 @@ Vec3 Elasticity::getBodyforce (const Vec3& X) const
 bool Elasticity::haveLoads () const
 {
   if (tracFld) return true;
+  if (fluxFld) return true;
   if (bodyFld) return true;
 
   for (unsigned short int i = 0; i < nsd; i++)
@@ -208,6 +212,7 @@ bool Elasticity::haveLoads () const
 
 void Elasticity::initIntegration (size_t, size_t nBp)
 {
+  tracVal.clear();
   tracVal.resize(nBp,std::make_pair(Vec3(),Vec3()));
 }
 

src/Integrands/Elasticity.h

@@ -45,6 +45,8 @@ public:
 
   //! \brief Defines the traction field to use in Neumann boundary conditions.
   void setTraction(TractionFunc* tf) { tracFld = tf; }
+  //! \brief Defines the traction field to use in Neumann boundary conditions.
+  void setTraction(VecFunc* tf) { fluxFld = tf; }
   //! \brief Defines the body force field.
   void setBodyForce(VecFunc* bf) { bodyFld = bf; }
 
@@ -216,7 +218,8 @@ protected:
   double    grav[3];  //!< Gravitation vector
 
   LocalSystem*  locSys;  //!< Local coordinate system for result output
-  TractionFunc* tracFld; //!< Pointer to boundary traction field
+  TractionFunc* tracFld; //!< Pointer to implicit boundary traction field
+  VecFunc*      fluxFld; //!< Pointer to explicit boundary traction field
   VecFunc*      bodyFld; //!< Pointer to body force field
 
   mutable std::vector<Vec3Pair> tracVal; //!< Traction field point values

src/Integrands/LinearElasticity.C

@@ -92,7 +92,7 @@ bool LinearElasticity::evalInt (LocalIntegral& elmInt, const FiniteElement& fe,
 bool LinearElasticity::evalBou (LocalIntegral& elmInt, const FiniteElement& fe,
 				const Vec3& X, const Vec3& normal) const
 {
-  if (!tracFld)
+  if (!tracFld && !fluxFld)
   {
     std::cerr <<" *** LinearElasticity::evalBou: No tractions."<< std::endl;
     return false;
@@ -107,11 +107,14 @@ bool LinearElasticity::evalBou (LocalIntegral& elmInt, const FiniteElement& fe,
   const double detJW = axiSymmetry ? 2.0*M_PI*X.x*fe.detJxW : fe.detJxW;
 
   // Evaluate the surface traction
-  Vec3 T = (*tracFld)(X,normal);
+  Vec3 T = this->getTraction(X,normal);
 
   // Store traction value for visualization
-  if (fe.iGP < tracVal.size())
+  if (fe.iGP < tracVal.size() && !T.isZero())
-    if (!T.isZero()) tracVal[fe.iGP] = std::make_pair(X,T);
+  {
+    tracVal[fe.iGP].first = X;
+    tracVal[fe.iGP].second += T;
+  }
 
   // Integrate the force vector
   Vector& ES = static_cast<ElmMats&>(elmInt).b[eS-1];

src/Integrands/Poisson.C

@@ -29,6 +29,7 @@ Poisson::Poisson (unsigned short int n) : nsd(n)
   kappa = 1.0;
 
   tracFld = 0;
+  fluxFld = 0;
   heatSrc = 0;
 
   // Only the current solution is needed
@@ -44,7 +45,12 @@ double Poisson::getHeat (const Vec3& X) const
 
 double Poisson::getTraction (const Vec3& X, const Vec3& n) const
 {
-  return tracFld ? (*tracFld)(X)*n : 0.0;
+  if (fluxFld)
+    return (*fluxFld)(X);
+  else if (tracFld)
+    return (*tracFld)(X)*n;
+  else
+    return 0.0;
 }
 
 
@@ -115,7 +121,7 @@ bool Poisson::evalInt (LocalIntegral& elmInt, const FiniteElement& fe,
 bool Poisson::evalBou (LocalIntegral& elmInt, const FiniteElement& fe,
 		       const Vec3& X, const Vec3& normal) const
 {
-  if (!tracFld)
+  if (!tracFld && !fluxFld)
   {
     std::cerr <<" *** Poisson::evalBou: No tractions."<< std::endl;
     return false;
@@ -128,16 +134,15 @@ bool Poisson::evalBou (LocalIntegral& elmInt, const FiniteElement& fe,
     return false;
   }
 
-  // Evaluate the heat flux q at point X
+  // Evaluate the Neumann value -q(X)*n
-  Vec3 q = (*tracFld)(X);
+  double trac = -this->getTraction(X,normal);
-
-  // Evaluate the Neumann value -q*n
-  double trac = -(q*normal);
 
   // Store traction value for visualization
-  if (fe.iGP < tracVal.size())
+  if (fe.iGP < tracVal.size() && abs(trac) > 1.0e8)
-    if (abs(trac) > 1.0e8)
+  {
-      tracVal[fe.iGP] = std::make_pair(X,trac*normal);
+    tracVal[fe.iGP].first = X;
+    tracVal[fe.iGP].second += trac*normal;
+  }
 
   // Integrate the Neumann value
   elMat.b.front().add(fe.N,trac*fe.detJxW);
@@ -148,6 +153,7 @@ bool Poisson::evalBou (LocalIntegral& elmInt, const FiniteElement& fe,
 
 void Poisson::initIntegration (size_t, size_t nBp)
 {
+  tracVal.clear();
   tracVal.resize(nBp,std::make_pair(Vec3(),Vec3()));
 }
 

src/Integrands/Poisson.h

@@ -35,6 +35,8 @@ public:
 
   //! \brief Defines the traction field to use in Neumann boundary conditions.
   void setTraction(VecFunc* tf) { tracFld = tf; }
+  //! \brief Defines the traction field to use in Neumann boundary conditions.
+  void setTraction(RealFunc* ff) { fluxFld = ff; }
   //! \brief Defines the heat source field.
   void setSource(RealFunc* src) { heatSrc = src; }
 
@@ -135,6 +137,7 @@ private:
 
 protected:
   VecFunc*  tracFld; //!< Pointer to boundary traction field
+  RealFunc* fluxFld; //!< Pointer to boundary normal flux field
   RealFunc* heatSrc; //!< Pointer to interior heat source
 
   mutable std::vector<Vec3Pair> tracVal; //!< Traction field point values

src/SIM/SIMbase.C

@@ -84,7 +84,7 @@ SIMbase::~SIMbase ()
     delete *i1;
 
   myModel.clear();
-  this->clearProperties();
+  this->SIMbase::clearProperties();
 
 #ifdef SP_DEBUG
   std::cout <<"Leaving SIMbase destructor"<< std::endl;
@@ -149,7 +149,7 @@ bool SIMbase::parseGeometryTag (const TiXmlElement* elem)
     int proc = 0;
     if (!utl::getAttribute(elem,"procs",proc))
       return false;
-    if (proc != nProc) // silently ignore
+    else if (proc != nProc) // silently ignore
       return true;
     if (myPid == 0)
       std::cout <<"\tNumber of partitions: "<< nProc << std::endl;
@@ -231,21 +231,36 @@ bool SIMbase::parseBCTag (const TiXmlElement* elem)
     }
   }
 
+  else if (!strcasecmp(elem->Value(),"neumann")) {
+    int code = 0, ndir = 0;
+    std::string type;
+    utl::getAttribute(elem,"code",code);
+    utl::getAttribute(elem,"direction",ndir);
+    utl::getAttribute(elem,"type",type,true);
+    if (elem->FirstChild()) {
+      std::cout <<"\tNeumann code "<< code <<" direction "<< ndir;
+      if (!type.empty()) std::cout <<" ("<< type <<")";
+      this->setNeumann(elem->FirstChild()->Value(),type,ndir,code);
+      std::cout << std::endl;
+    }
+  }
+
   else if (!strcasecmp(elem->Value(),"dirichlet") && !ignoreDirichlet) {
     int code = 0;
+    std::string type;
     utl::getAttribute(elem,"code",code);
-    double val = elem->FirstChild() ? atof(elem->FirstChild()->Value()) : 0.0;
+    utl::getAttribute(elem,"type",type,true);
-    if (val == 0.0)
+    const TiXmlNode* dval = elem->FirstChild();
+    if (!dval || (atof(dval->Value()) == 0.0 && type != "expression")) {
       setPropertyType(code,Property::DIRICHLET);
-    else {
+      std::cout <<"\tDirichlet code "<< code <<": (fixed)"<< std::endl;
+    }
+    else if (dval) {
       setPropertyType(code,Property::DIRICHLET_INHOM);
-      std::string function;
+      std::cout <<"\tDirichlet code "<< code;
-      if (utl::getAttribute(elem,"function",function)) {
+      if (!type.empty()) std::cout <<" ("<< type <<")";
-        char* func = strtok(const_cast<char*>(function.c_str())," ");
+      myScalars[code] = utl::parseRealFunc(dval->Value(),type);
-        myScalars[code] = const_cast<RealFunc*>(utl::parseRealFunc(func,val));
+      std::cout << std::endl;
-      }
-      else
-        myScalars[code] = new ConstFunc(val);
     }
   }
 
@@ -261,7 +276,7 @@ bool SIMbase::parseOutputTag (const TiXmlElement* elem)
     return opt.parseOutputTag(elem);
 
   const TiXmlElement* point = elem->FirstChildElement("point");
-  for (int i = 1; point; i++)
+  for (int i = 1; point; i++, point = point->NextSiblingElement())
   {
     int patch = 0;
     ResultPoint thePoint;
@@ -276,7 +291,6 @@ bool SIMbase::parseOutputTag (const TiXmlElement* elem)
     if (utl::getAttribute(point,"w",thePoint.par[2]) && myPid == 0)
       std::cout <<' '<< thePoint.par[2];
     myPoints.push_back(thePoint);
-    point = point->NextSiblingElement();
   }
   if (myPid == 0)
     std::cout << std::endl;
@@ -824,6 +838,38 @@ bool SIMbase::setVecProperty (int code, Property::Type ptype, VecFunc* field)
 }
 
 
+bool SIMbase::setNeumann (const std::string& prop, const std::string& type,
+			  int direction, int code)
+{
+  if (direction == 0 && this->getNoFields() == 1)
+  {
+    RealFunc* f = utl::parseRealFunc(prop,type);
+    if (f)
+      myScalars[code] = f;
+    else
+      return false;
+  }
+  else if (direction < 0 || this->getNoFields() == 1)
+  {
+    VecFunc* f = utl::parseVecFunc(prop,type);
+    if (f)
+      myVectors[code] = f;
+    else
+      return false;
+  }
+  else
+  {
+    TractionFunc* f = utl::parseTracFunc(prop,type,direction);
+    if (f)
+      myTracs[code] = f;
+    else
+      return false;
+  }
+
+  return this->setPropertyType(code,Property::NEUMANN);
+}
+
+
 bool SIMbase::initSystem (int mType, size_t nMats, size_t nVec)
 {
   if (!mySam) return false;
@@ -2154,7 +2200,7 @@ bool SIMbase::dumpResults (const Vector& psol, double time, std::ostream& os,
     for (j = 0; j < points.size(); j++)
     {
       if (!formatted)
-        os << time <<" ";
+	os << time <<" ";
       else if (points[j] < 0)
 	os <<"  Point #"<< -points[j] <<":\tsol1 =";
       else
@@ -2168,8 +2214,8 @@ bool SIMbase::dumpResults (const Vector& psol, double time, std::ostream& os,
 
       if (opt.discretization >= ASM::Spline)
       {
-        if (formatted && sol2.rows() > 0)
+	if (formatted && sol2.rows() > 0)
-          os <<"\n\t\tsol2 =";
+	  os <<"\n\t\tsol2 =";
 	for (k = 1; k <= sol2.rows(); k++)
 	  os << std::setw(flWidth) << utl::trunc(sol2(k,j+1));
       }
@@ -2178,8 +2224,8 @@ bool SIMbase::dumpResults (const Vector& psol, double time, std::ostream& os,
 	// Print nodal reaction forces for nodes with prescribed DOFs
 	if (mySam->getNodalReactions(points[j],*reactionForces,reactionFS))
 	{
-          if (formatted)
+	  if (formatted)
-            os <<"\n\t\treac =";
+	    os <<"\n\t\treac =";
 	  for (k = 0; k < reactionFS.size(); k++)
 	    os << std::setw(flWidth) << utl::trunc(reactionFS[k]);
 	}
@@ -2221,35 +2267,35 @@ bool SIMbase::project (Matrix& ssol, const Vector& psol,
     switch (pMethod) {
     case SIMoptions::GLOBAL:
       if (msgLevel > 1 && i == 0)
-	std::cout <<"\tGreville point projection"<< std::endl;
+        std::cout <<"\tGreville point projection"<< std::endl;
       if (!myModel[i]->evalSolution(values,*myProblem))
-	return false;
+        return false;
       break;
 
     case SIMoptions::DGL2:
       if (msgLevel > 1 && i == 0)
-	std::cout <<"\tDiscrete global L2-projection"<< std::endl;
+        std::cout <<"\tDiscrete global L2-projection"<< std::endl;
       if (!myModel[i]->globalL2projection(values,*myProblem))
-	return false;
+        return false;
       break;
 
     case SIMoptions::CGL2:
       if (msgLevel > 1 && i == 0)
-	std::cout <<"\tContinuous global L2-projection"<< std::endl;
+        std::cout <<"\tContinuous global L2-projection"<< std::endl;
       if (!myModel[i]->globalL2projection(values,*myProblem,true))
-	return false;
+        return false;
       break;
 
     case SIMoptions::SCR:
       if (msgLevel > 1 && i == 0)
-	std::cout <<"\tSuperconvergent recovery"<< std::endl;
+        std::cout <<"\tSuperconvergent recovery"<< std::endl;
       if (!myModel[i]->evalSolution(values,*myProblem,NULL,'S'))
-	return false;
+        return false;
       break;
 
     case SIMoptions::VDSA:
       if (msgLevel > 1 && i == 0)
-        std::cout <<"\tVDSA projection"<< std::endl;
+        std::cout <<"\tVariation diminishing projection"<< std::endl;
       if (!myModel[i]->evalSolution(values,*myProblem,NULL,'A'))
         return false;
       break;
@@ -2270,7 +2316,7 @@ bool SIMbase::project (Matrix& ssol, const Vector& psol,
 
     default:
       std::cerr <<" *** SIMbase::project: Projection method "<< pMethod
-		<<" not implemented."<< std::endl;
+                <<" not implemented."<< std::endl;
       return false;
     }
 

src/SIM/SIMbase.h

@@ -488,6 +488,14 @@ protected:
   //! \param[in] pindex 0-based index into problem-dependent property container
   bool setPropertyType(int code, Property::Type ptype, int pindex = -1);
 
+  //! \brief Defines a Neumann boundary condition property by parsing a string.
+  //! \param[in] prop The string to parse the property definition from
+  //! \param[in] type Additional option defining the type of property definition
+  //! \param[in] ndir Direction of the surface traction on the Neumann boundary
+  //! \param[in] code The property code to be associated with this property
+  bool setNeumann(const std::string& prop, const std::string& type,
+		  int ndir, int code);
+
   //! \brief Preprocesses a user-defined Dirichlet boundary property.
   //! \param[in] patch 1-based index of the patch to receive the property
   //! \param[in] lndx Local index of the boundary item to receive the property