Added: Additional MADOF support in SIMbase for mixed problems.

This allows DOF vectors with an "unusual" number of DOFs on a given basis.
Added usage of this in the HDF5Writer for such fields.

src/SIM/SIMbase.C

@@ -1488,7 +1488,7 @@ size_t SIMbase::getNoRHS () const
 }
 
 
-char SIMbase::getNoBasis () const
+unsigned char SIMbase::getNoBasis () const
 {
   size_t result = myModel.empty() ? 0 : myModel.front()->getNoBasis();
 #ifdef SP_DEBUG
@@ -2509,23 +2509,39 @@ bool SIMbase::extractPatchSolution (IntegrandBase* problem,
 
 
 size_t SIMbase::extractPatchSolution (const Vector& sol, Vector& vec,
-                                      int pindx, unsigned char nndof) const
+                                      int pindx, unsigned char nndof,
+                                      unsigned char basis) const
 {
   ASMbase* pch = pindx >= 0 ? this->getPatch(pindx+1) : nullptr;
   if (!pch || sol.empty()) return 0;
 
-  pch->extractNodeVec(sol,vec,nndof);
+  if (basis != 0 && nndof != 0 &&
+      nndof != this->getNoFields(basis) && this->getNoFields(2) > 0)
+  {
+    // Need an additional MADOF
+    const std::vector<int>& madof = this->getMADOF(basis,nndof);
+    pch->extractNodalVec(sol,vec,madof.data(),madof.size());
+  }
+  else
+    pch->extractNodeVec(sol,vec,nndof,basis);
 
-  return (nndof > 0 ? nndof : pch->getNoFields(1))*pch->getNoNodes(1);
+  return vec.size();
 }
 
 
 bool SIMbase::injectPatchSolution (Vector& sol, const Vector& vec,
-                                   int pindx, unsigned char nndof) const
+                                   int pindx, unsigned char nndof,
+                                   unsigned char basis) const
 {
   ASMbase* pch = pindx >= 0 ? this->getPatch(pindx+1) : nullptr;
+  if (!pch) return false;
 
-  return pch ? pch->injectNodeVec(vec,sol,nndof) : false;
+  if (basis > 0 && nndof > 0 &&
+      nndof != this->getNoFields(basis) && this->getNoFields(2) > 0)
+    // Need an additional MADOF
+    return pch->injectNodalVec(vec,sol,this->getMADOF(basis,nndof),basis);
+  else
+    return pch->injectNodeVec(vec,sol,nndof);
 }
 
 
@@ -2614,3 +2630,45 @@ bool SIMbase::refine (const LR::RefineData& prm,
 
   return isRefined;
 }
+
+
+bool SIMbase::addMADOF (unsigned char basis, unsigned char nndof)
+{
+  int key = basis << 16 + nndof;
+  auto it = mixedMADOFs.find(key);
+  if (it != mixedMADOFs.end())
+    return false; // This MADOF already calculated
+
+  std::vector<int>& madof = mixedMADOFs[key];
+  madof.resize(this->getNoNodes(true)+1,0);
+
+  char nType = basis <= 1 ? 'D' : 'P' + basis-2;
+  for (size_t i = 0; i < myModel.size(); i++)
+    for (size_t j = 0; j < myModel[i]->getNoNodes(); j++)
+    {
+      int n = myModel[i]->getNodeID(j+1);
+      if (n > 0 && myModel[i]->getNodeType(j+1) == nType)
+        madof[n] = nndof;
+      else if (n > 0)
+        madof[n] = myModel[i]->getNodalDOFs(j+1);
+    }
+
+  madof[0] = 1;
+  for (size_t n = 1; n < madof.size(); n++)
+    madof[n] += madof[n-1];
+
+  return true;
+}
+
+
+const std::vector<int>& SIMbase::getMADOF (unsigned char basis,
+                                           unsigned char nndof) const
+{
+  int key = basis << 16 + nndof;
+  auto it = mixedMADOFs.find(key);
+  if (it != mixedMADOFs.end())
+    return it->second;
+
+  static std::vector<int> empty;
+  return empty;
+}

src/SIM/SIMbase.h

@@ -228,7 +228,7 @@ public:
   //! \brief Returns the number of right-hand-side vectors.
   virtual size_t getNoRHS() const;
   //! \brief Returns the number of bases in the model.
-  char getNoBasis() const;
+  unsigned char getNoBasis() const;
 
   //! \brief Returns the type (DOF classification) of the specified global node.
   char getNodeType(int inod) const;
@@ -588,17 +588,21 @@ public:
   //! \param[out] vec Local solution vector associated with specified patch
   //! \param[in] pindx Local patch index to extract solution vector for
   //! \param[in] nndof Number of DOFs per node (optional)
+  //! \param[in] basis Basis to extract for (optional)
   //! \return Total number of DOFs in the patch (first basis only if mixed)
   size_t extractPatchSolution(const Vector& sol, Vector& vec, int pindx,
-                              unsigned char nndof = 0) const;
+                              unsigned char nndof = 0,
+                              unsigned char basis = 0) const;
 
   //! \brief Injects a patch-wise solution vector into the global vector.
   //! \param sol Global primary solution vector in DOF-order
   //! \param[in] vec Local solution vector associated with specified patch
   //! \param[in] pindx Local patch index to inject solution vector for
   //! \param[in] nndof Number of DOFs per node (optional)
+  //! \param[in] basis Basis to inject for (optional)
   bool injectPatchSolution(Vector& sol, const Vector& vec, int pindx,
-                           unsigned char nndof = 0) const;
+                           unsigned char nndof = 0,
+                           unsigned char basis = 0) const;
 
   //! \brief Extracts element results for a specified patch.
   //! \param[in] glbRes Global element result array
@@ -694,6 +698,18 @@ protected:
   //! \param[in] silence If \e true, suppress threading group outprint
   void generateThreadGroups(const Property& p, bool silence = false);
 
+  //! \brief Adds a MADOF with an extraordinary number of DOFs on a given basis.
+  //! \param[in] basis The basis to specify number of DOFs for
+  //! \param[in] nndof Number of nodal DOFs on the given basis
+  bool addMADOF(unsigned char basis, unsigned char nndof);
+
+private:
+  //! \brief Returns an extraordinary MADOF array.
+  //! \param[in] basis The basis to specify number of DOFs for
+  //! \param[in] nndof Number of nodal DOFs on the given basis
+  const std::vector<int>& getMADOF(unsigned char basis,
+                                   unsigned char nndof) const;
+
 public:
   static bool ignoreDirichlet; //!< Set to \e true for free vibration analysis
   static bool preserveNOrder;  //!< Set to \e true to preserve node ordering
@@ -751,6 +767,9 @@ protected:
 private:
   size_t nIntGP; //!< Number of interior integration points in the whole model
   size_t nBouGP; //!< Number of boundary integration points in the whole model
+
+  //! Additional MADOF arrays for mixed problems (extraordinary DOF counts)
+  std::map<int, std::vector<int> > mixedMADOFs;
 };
 
 #endif

src/SIM/Test/TestSIM.C

@@ -26,6 +26,11 @@ public:
   {
     Dim::myProblem = new TestProjectIntegrand(Dim::dimension);
   }
+
+  bool addMixedMADOF(unsigned char basis, unsigned char nndof)
+  {
+    return this->addMADOF(basis, nndof);
+  }
 private:
   class TestProjectIntegrand : public IntegrandBase {
   public:
@@ -129,3 +134,42 @@ TEST(TestSIM3D, ProjectSolutionMixed)
       for (size_t i = 0; i < 3; ++i)
         ASSERT_FLOAT_EQ(ssol(1, n++), i/2.0 + j/2.0 + k/2.0);
 }
+
+
+TEST(TestSIM, InjectPatchSolution)
+{
+  TestProjectSIM<SIM2D> sim({1,1});
+  sim.createDefaultModel();
+  ASSERT_TRUE(sim.preprocess());
+
+  Vector sol(2*sim.getNoNodes(true, 1) + sim.getNoNodes(true, 2));
+  Vector lsol(2*sim.getNoNodes(true, 1));
+  for (size_t i = 0; i < sim.getNoNodes(true,1); ++i)
+    lsol[2*i] = lsol[2*i+1] = i+1;
+
+  ASSERT_TRUE(sim.addMixedMADOF(1, 2));
+  sim.injectPatchSolution(sol, lsol, 0, 2, 1);
+  size_t ofs = 0;
+  for (size_t i = 0; i < sim.getNoNodes(true,1); ++i, ofs += 2) {
+    ASSERT_FLOAT_EQ(sol[ofs], i+1);
+    ASSERT_FLOAT_EQ(sol[ofs+1], i+1);
+  }
+  for (size_t i = 0; i < sim.getNoNodes(true,2); ++i, ++ofs)
+    ASSERT_FLOAT_EQ(sol[ofs], 0);
+
+  ASSERT_TRUE(sim.addMixedMADOF(2, 2));
+  Vector sol2(sim.getNoNodes(true, 1) + 2*sim.getNoNodes(true, 2));
+  Vector lsol2(2*sim.getNoNodes(true, 2));
+  for (size_t i = 0; i < sim.getNoNodes(true,2); ++i)
+    lsol2[2*i] = lsol2[2*i+1] = i+1;
+
+  sim.injectPatchSolution(sol2, lsol2, 0, 2, 2);
+  ofs = 0;
+  for (size_t i = 0; i < sim.getNoNodes(true,1); ++i, ++ofs)
+    ASSERT_FLOAT_EQ(sol2[ofs], 0);
+
+  for (size_t i = 0; i < sim.getNoNodes(true,2); ++i, ofs += 2) {
+    ASSERT_FLOAT_EQ(sol2[ofs], i+1);
+    ASSERT_FLOAT_EQ(sol2[ofs+1], i+1);
+  }
+}

src/Utility/HDF5Writer.C

@@ -397,7 +397,7 @@ void HDF5Writer::writeSIM (int level, const DataEntry& entry,
   if (level == 0 || geometryUpdated || (abs(entry.second.results) & DataExporter::GRID)) {
     writeBasis(sim,basisname,1,level);
     if (sim->mixedProblem())
-      for (size_t b=2; b <= sim->getPatch(1)->getNoBasis(); ++b) {
+      for (size_t b=2; b <= sim->getNoBasis(); ++b) {
         std::stringstream str;
         str << sim->getName() << "-" << b;
         writeBasis(sim,str.str(),b,level);
@@ -445,21 +445,27 @@ void HDF5Writer::writeSIM (int level, const DataEntry& entry,
       if (abs(results) & DataExporter::PRIMARY) {
         Vector psol;
         int ncmps = entry.second.ncmps;
-        size_t ndof1 = sim->extractPatchSolution(*sol,psol,loc-1,ncmps);
+        if (entry.second.results < 0) { // field assumed to be on basis 1 for now
-        if (sim->mixedProblem())
+          size_t ndof1 = sim->extractPatchSolution(*sol, psol, loc-1, ncmps, 1);
-        {
+          writeArray(group2, entry.second.description,
-          size_t ofs = 0;
+                     ndof1, psol.ptr(), H5T_NATIVE_DOUBLE);
-          for (size_t b=1; b <= sim->getPatch(loc)->getNoBasis(); ++b) {
+        } else {
-            ndof1 = sim->getPatch(loc)->getNoNodes(b)*sim->getPatch(loc)->getNoFields(b);
+          size_t ndof1 = sim->extractPatchSolution(*sol,psol,loc-1,ncmps);
-            writeArray(group2,prefix+prob->getField1Name(10+b),ndof1,
+          if (sim->mixedProblem())
-                       psol.ptr()+ofs,H5T_NATIVE_DOUBLE);
+          {
-            ofs += ndof1;
+            size_t ofs = 0;
+            for (size_t b=1; b <= sim->getNoBasis(); ++b) {
+              ndof1 = sim->getPatch(loc)->getNoNodes(b)*sim->getPatch(loc)->getNoFields(b);
+              writeArray(group2,prefix+prob->getField1Name(10+b),ndof1,
+                         psol.ptr()+ofs,H5T_NATIVE_DOUBLE);
+              ofs += ndof1;
+            }
+          }
+          else {
+            writeArray(group2, usedescription ? entry.second.description:
+                                                prefix+prob->getField1Name(11),
+                                         ndof1, psol.ptr(), H5T_NATIVE_DOUBLE);
           }
-        }
-        else {
-          writeArray(group2, usedescription ? entry.second.description:
-                                              prefix+prob->getField1Name(11),
-                                       ndof1, psol.ptr(), H5T_NATIVE_DOUBLE);
         }
       }
 
@@ -533,10 +539,13 @@ void HDF5Writer::writeSIM (int level, const DataEntry& entry,
                    0, &dummy, H5T_NATIVE_DOUBLE);
       }
       if (abs(results) & DataExporter::PRIMARY) {
-        if (sim->mixedProblem())
+        if (entry.second.results < 0) {
+          writeArray(group2, entry.second.description,
+                     0, &dummy, H5T_NATIVE_DOUBLE);
+        }
+        else if (sim->mixedProblem())
         {
-          writeArray(group2,prefix+entry.first,0,&dummy,H5T_NATIVE_DOUBLE);
+          for (size_t b=1; b <= sim->getNoBasis(); ++b)
-          for (size_t b=1; b <= sim->getPatch(loc)->getNoBasis(); ++b)
             writeArray(group2,prefix+prob->getField1Name(10+b),0,&dummy,H5T_NATIVE_DOUBLE);
         }
         else

src/Utility/XMLWriter.C

@@ -222,7 +222,10 @@ void XMLWriter::writeSIM (int level, const DataEntry& entry, bool,
   }
 
   if (results & DataExporter::PRIMARY) {
-    if (sim->mixedProblem())
+    if (entry.second.results < 0)
+      addField(entry.second.description, entry.second.description,
+               basisname, cmps, sim->getNoPatches(), "field");
+    else if (sim->mixedProblem())
     {
       // primary solution vector
       addField(prefix+entry.first,entry.second.description,basisname,