Fixed: Missing argument to NextSiblingElement() calls.

Added: Evaluation of 2ndary solution at specified Cartesian points.

Author: kmokstad

src/SIM/SIM1D.C

@@ -148,7 +148,7 @@ bool SIM1D::parseGeometryTag (const tinyxml2::XMLElement* elem)
     utl::getAttribute(elem,"offset",offset);
 
     const tinyxml2::XMLElement* child = elem->FirstChildElement("connection");
-    for (; child; child = child->NextSiblingElement())
+    for (; child; child = child->NextSiblingElement("connection"))
     {
       ASM::Interface ifc;
       if (utl::getAttribute(child,"master",ifc.master))

src/SIM/SIM2D.C

@@ -187,7 +187,7 @@ bool SIM2D::parseGeometryTag (const tinyxml2::XMLElement* elem)
 
     std::vector<Interface> top;
     const tinyxml2::XMLElement* child = elem->FirstChildElement("connection");
-    for (; child; child = child->NextSiblingElement())
+    for (; child; child = child->NextSiblingElement("connection"))
     {
       ASM::Interface ifc;
       bool rever = false, periodic = false;

src/SIM/SIM3D.C

@@ -178,7 +178,7 @@ bool SIM3D::parseGeometryTag (const tinyxml2::XMLElement* elem)
     utl::getAttribute(elem,"offset",offset);
 
     const tinyxml2::XMLElement* child = elem->FirstChildElement("connection");
-    for (; child; child = child->NextSiblingElement())
+    for (; child; child = child->NextSiblingElement("connection"))
     {
       ASM::Interface ifc;
       bool periodic = false;

src/SIM/SIMinput.C

@@ -451,7 +451,7 @@ bool SIMinput::parseBCTag (const tinyxml2::XMLElement* elem)
   else if (!strcasecmp(elem->Value(),"propertycodes"))
   {
     const tinyxml2::XMLElement* code = elem->FirstChildElement("code");
-    for (; code; code = code->NextSiblingElement())
+    for (; code; code = code->NextSiblingElement("code"))
     {
       int icode = 0;
       utl::getAttribute(code,"value",icode);

src/SIM/SIMoutput.C

@@ -106,8 +106,7 @@ bool SIMoutput::parseOutputTag (const tinyxml2::XMLElement* elem)
 {
   IFEM::cout <<"  Parsing <"<< elem->Value() <<">"<< std::endl;
 
-  const char* funcval = utl::getValue(elem,"function");
-  if (funcval)
+  if (const char* funcval = utl::getValue(elem,"function"); funcval)
   {
     std::string name;
     utl::getAttribute(elem,"name",name);
@@ -141,7 +140,7 @@ bool SIMoutput::parseOutputTag (const tinyxml2::XMLElement* elem)
   // Parse the result point specifications.
   // Can either be explicit points, lines or a grid.
   const tinyxml2::XMLElement* point = elem->FirstChildElement("point");
-  for (int i = 1; point; i++, point = point->NextSiblingElement())
+  for (int i = 1; point; i++, point = point->NextSiblingElement("point"))
   {
     ResultPoint thePoint(3);
     if (int patch = 0; utl::getAttribute(point,"patch",patch) && patch > 0)
@@ -150,21 +149,30 @@ bool SIMoutput::parseOutputTag (const tinyxml2::XMLElement* elem)
     IFEM::cout <<"\tPoint "<< i <<": P"<< thePoint.patch;
     if (utl::getAttribute(point,"node",thePoint.inod) && thePoint.inod > 0)
       IFEM::cout <<" node = "<< thePoint.inod;
+    else if (point->FirstChild() &&
+             utl::parseVec(thePoint.X,point->FirstChild()->Value()))
+    {
+      // Cartesian coordinates of a spatial point is specified,
+      thePoint.iel = -1; // need to search for the matching element/node
+      IFEM::cout <<" X = "<< thePoint.X;
+    }
     else
+    {
       IFEM::cout <<" xi =";
-    if (utl::getAttribute(point,"u",thePoint.u[0]))
-      IFEM::cout <<' '<< thePoint.u[0];
-    if (utl::getAttribute(point,"v",thePoint.u[1]))
-      IFEM::cout <<' '<< thePoint.u[1];
-    if (utl::getAttribute(point,"w",thePoint.u[2]))
-      IFEM::cout <<' '<< thePoint.u[2];
+      if (utl::getAttribute(point,"u",thePoint.u[0]))
+        IFEM::cout <<' '<< thePoint.u[0];
+      if (utl::getAttribute(point,"v",thePoint.u[1]))
+        IFEM::cout <<' '<< thePoint.u[1];
+      if (utl::getAttribute(point,"w",thePoint.u[2]))
+        IFEM::cout <<' '<< thePoint.u[2];
+    }
     IFEM::cout << std::endl;
 
     addPoint(thePoint);
   }
 
   const tinyxml2::XMLElement* line = elem->FirstChildElement("line");
-  for (int j = 1; line; j++, line = line->NextSiblingElement())
+  for (int j = 1; line; j++, line = line->NextSiblingElement("line"))
   {
     ResultPoint thePoint(3);
     if (int patch = 0; utl::getAttribute(line,"patch",patch) && patch > 0)
@@ -206,26 +214,22 @@ bool SIMoutput::parseOutputTag (const tinyxml2::XMLElement* elem)
   }
 
   for (const tinyxml2::XMLElement* child = elem->FirstChildElement("elements");
-       child; child = child->NextSiblingElement())
+       child; child = child->NextSiblingElement("elements"))
     if (const tinyxml2::XMLNode* elist = child->FirstChild(); elist)
     {
       ResultPoint thePoint;
       if (int patch = 0; utl::getAttribute(child,"patch",patch) && patch > 0)
         thePoint.patch = patch;
 
-      IntVec elements;
-      utl::parseIntegers(elements,elist->Value());
-      for (int iel : elements)
-      {
-        thePoint.iel = iel;
-        addPoint(thePoint);
-      }
-
-      if (!elements.empty())
+      if (IntVec elements; utl::parseIntegers(elements,elist->Value()))
       {
         IFEM::cout <<"\tElement center output: P"<< thePoint.patch <<" (";
-        for (const ResultPoint& rp : myPoints.back().second)
-          IFEM::cout <<" "<< rp.iel;
+        for (int iel : elements)
+        {
+          IFEM::cout <<" "<< iel;
+          thePoint.iel = iel;
+          addPoint(thePoint);
+        }
         IFEM::cout <<" )"<< std::endl;
       }
     }
@@ -236,7 +240,7 @@ bool SIMoutput::parseOutputTag (const tinyxml2::XMLElement* elem)
   else if (grid)
     idxGrid = myPoints.size();
 
-  for (int g = 1; grid; g++, grid = grid->NextSiblingElement())
+  for (int g = 1; grid; g++, grid = grid->NextSiblingElement("grid"))
   {
     ResultPoint thePoint;
     if (int patch = 0; utl::getAttribute(grid,"patch",patch) && patch > 0)
@@ -413,7 +417,7 @@ void SIMoutput::preprocessResPtGroup (std::string& ptFile, ResPointVec& points)
 {
   if (points.empty()) return;
 
-  // Check if we have Cartesian grid input
+  // Check if we have a regular Cartesian grid input
   size_t nX = 0;
   IntMat pointMap;
   if (points.front().inod < 0 && points.back().inod < 0)
@@ -426,23 +430,52 @@ void SIMoutput::preprocessResPtGroup (std::string& ptFile, ResPointVec& points)
       nX = 0;
   }
 
+  // Check if we have general spatial points input (Cartesian coordinates)
+  using PatchPoints = std::pair<std::vector<ASMbase::PointParams>,Vec3Vec>;
+  std::map<ASMbase*,PatchPoints> pchPoints;
+  if (pointMap.empty())
+    for (ResultPoint& pt : points)
+      if (pt.iel < 0 || pt.inod < 0)
+        if (ASMbase* pch = this->getPatch(pt.patch,true); pch && !pch->empty())
+        {
+          pchPoints[pch].second.emplace_back(pt.X);
+          pt.iel = -pchPoints[pch].second.size();
+        }
+
+  // Find element and/or local parameters of the spatial points
+  for (std::pair<ASMbase* const,PatchPoints>& pch : pchPoints)
+    if (!pch.first->findPoints(pch.second.second,pch.second.first))
+      IFEM::cout <<"  ** Failed to locate some result points."<< std::endl;
+
   int iPoint = 0;
   size_t iX = 0, iY = 0;
   for (ResPointVec::iterator pit = points.begin(); pit != points.end();)
   {
     bool pointIsOK = false;
-    ASMbase* pch = this->getPatch(pit->patch,true);
-    if (pch && !pch->empty())
+    if (ASMbase* pch = this->getPatch(pit->patch,true); pch && !pch->empty())
     {
       if ((pointIsOK = pit->inod > 0)) // A nodal number is specified
         pit->X = pch->getCoord(pit->inod);
 
-      else if (pit->inod < 0) // A spatial point is specified
-        pointIsOK = fabs(pch->findPoint(pit->X,pit->u)) < 1.0e-3;
-
+      else if (pit->iel < 0 && !pchPoints.empty())
+      {
+        // A spatial point is specified, use the point searching results
+        size_t iP = -pit->iel-1;
+        std::map<ASMbase*,PatchPoints>::const_iterator it = pchPoints.find(pch);
+        if (it != pchPoints.end() && iP < it->second.first.size())
+          if ((pointIsOK = fabs(it->second.first[iP].dist) < 1.0e-3))
+          {
+            pit->npar = pch->getNoParamDim();
+            pit->X = it->second.second[iP];
+            pit->iel = it->second.first[iP].iel;
+            memcpy(pit->u,it->second.first[iP].u,sizeof(double)*3);
+          }
+      }
       else if (pit->npar > 0) // A parametric point is specified
+      {
+        pit->npar = pch->getNoParamDim();
         pointIsOK = (pit->inod = pch->evalPoint(pit->u,pit->u,pit->X)) >= 0;
-
+      }
       else if (pit->iel < 1) // All nodes or element centers
         pointIsOK = true;
 
@@ -455,8 +488,6 @@ void SIMoutput::preprocessResPtGroup (std::string& ptFile, ResPointVec& points)
 
     if (pointIsOK)
     {
-      if (pit->npar > 0)
-        pit->npar = pch->getNoParamDim();
       ++pit;
       ++iPoint;
       if (nX > 0)
@@ -498,17 +529,33 @@ void SIMoutput::preprocessResPtGroup (std::string& ptFile, ResPointVec& points)
         IFEM::cout <<", node #"<< pt.inod;
       else if (pt.iel > 0)
         IFEM::cout <<", element #"<< pt.iel;
-      if (pt.inod > 0 && myModel.size() > 1)
+      ASMbase* pch = this->getPatch(pt.patch,true);
+      if (pch)
       {
-        ASMbase* pch = this->getPatch(pt.patch,true);
-        IFEM::cout <<", global #"<< pch->getNodeID(pt.inod);
+        int globalId = 0;
+        if (pt.inod > 0)
+        {
+          if (int nodeId = pch->getNodeID(pt.inod);
+              nodeId != pt.inod || myModel.size() > 1)
+            globalId = nodeId;
+        }
+        else if (pt.iel > 0)
+        {
+          if (int elmId = pch->getElmID(pt.iel);
+              elmId != pt.iel || myModel.size() > 1)
+            globalId = elmId;
+        }
+        if (globalId > 0)
+          IFEM::cout <<", global #"<< globalId;
       }
       if (pt.npar > 0 || pt.iel > 0)
         IFEM::cout <<", X = "<< pt.X << std::endl;
       else if (pt.npar < 0)
         IFEM::cout <<", element centers will be used."<< std::endl;
       else
         IFEM::cout <<", nodal points will be used."<< std::endl;
+      if (pch && pt.iel > 0)
+        pch->printElmInfo(pt.iel,myProblem);
     }
   }
 
@@ -539,7 +586,7 @@ void SIMoutput::preprocessResPtGroup (std::string& ptFile, ResPointVec& points)
   for (const ResultPoint& pt : points)
   {
     os << 0.0 <<" "; // dummy time
-    for (short int k = 0; k < pt.npar; k++)
+    for (short int k = 0; k < nsd; k++)
       os << std::setw(flWidth) << pt.X[k];
     os << std::endl;
   }
@@ -561,8 +608,7 @@ bool SIMoutput::merge (SIMbase* that, const std::map<int,int>* old2new,
     for (ResultPoint& rp : rptp.second)
     {
       if (++iPoint == 1) IFEM::cout <<'\n';
-      ASMbase* pch = that->getPatch(rp.patch,true);
-      if (pch && rp.inod > 0)
+      if (ASMbase* pch = that->getPatch(rp.patch,true); pch && rp.inod > 0)
         IFEM::cout <<"Result point #"<< iPoint <<": patch #"<< rp.patch
                    <<",  local node "<< rp.inod
                    <<"  global node "<< pch->getNodeID(rp.inod) << std::endl;
@@ -1796,11 +1842,11 @@ bool SIMoutput::dumpMatlabGrid (std::ostream& os, const std::string& name,
   os <<"]="<< name << std::endl;
 
   // Write out nodes for the specified topology sets
-  TopologySet::const_iterator tit;
   for (const std::string& setname : sets)
   {
     std::set<int> nodeSet;
-    if ((tit = myEntitys.find(setname)) != myEntitys.end())
+    if (TopologySet::const_iterator tit = myEntitys.find(setname);
+        tit != myEntitys.end())
       for (const TopItem& top : tit->second)
         if (ASMbase* pch = this->getPatch(top.patch); pch)
           if (top.idim+1 == pch->getNoParamDim())
@@ -2021,7 +2067,7 @@ bool SIMoutput::dumpResults (const Vector& psol, double time,
       }
 
       bool newLine = false;
-      // Convenience lambda returning proper newline character before identifier.
+      // Convenience lambda returning proper newline char before identifier.
       auto&& newLineChr = [&newLine]() { return newLine ? "\n\t\t" : ":\t"; };
 
       if (formatted && !sol1.empty())
@@ -2124,6 +2170,8 @@ bool SIMoutput::evalResults (const Vectors& psol, const ResPointVec& gPoints,
       else if (pt.iel > 0)
       {
         elms.push_back(pt.iel);
+        for (short int k = 0; k < pt.npar; k++)
+          params[k].push_back(pt.u[k]);
       }
       else if (pt.npar < 0)
       {
@@ -2226,7 +2274,8 @@ bool SIMoutput::evalResults (const Vectors& psol, const ResPointVec& gPoints,
   if (elms.empty())
     ok = patch->evalSolution(tmp,*myProblem,params.data(),points.empty());
   else // evaluate at specified elements
-    ok = patch->evalSolution(tmp,*myProblem,elms);
+    ok = patch->evalSolution(tmp,*myProblem,elms,
+                             params.front().empty() ? nullptr : params.data());
 
   if ((ok &= augment(sol2,tmp)) && getNames)
     for (size_t i = 0; i < myProblem->getNoFields(2); i++)

src/SIM/TimeStep.C

@@ -118,7 +118,7 @@ bool TimeStep::parse (const tinyxml2::XMLElement* elem)
   if (f2 > 1.0) f2 = 1.0;
 
   const tinyxml2::XMLElement* child = elem->FirstChildElement("step");
-  for (; child; child = child->NextSiblingElement())
+  for (; child; child = child->NextSiblingElement("step"))
   {
     double start = 0.0, end = 0.0;
     utl::getAttribute(child,"start",start);