Subdivision for mixed bases

Apps/Common/MultiPatchModelGenerator.h

@@ -20,6 +20,23 @@
 #include <GoTools/trivariate/SplineVolume.h>
 
 
+//! \brief Extend basis (raise order and continuity) tuned for subdivided patches.
+//! \param[in] old Spline basis to extend
+//! \return New spline basis
+Go::BsplineBasis extendedBasis(const Go::BsplineBasis& old);
+//! \brief Helper for parsing periodic_? variable
+//! \param[in] nb Number of bases
+//! \param[in] basis Basis index under scrutiny
+//! \param[in] periodic Periodicity information to parse
+//! \return true if basis is periodic, false otherwise
+//! \return New spline basis
+static bool basisIsIn(size_t nb, int basis, int periodic) {
+  for (size_t p=0; p<nb; ++p)
+    if (basis == (periodic/(int)pow(10,p))%10)
+      return true;
+  return false; };
+
+
 /*!
   \brief 1D multi-patch model generator for FEM simulators.
   \details Generate a line split in a given number of blocks.
@@ -94,6 +111,9 @@ class MultiPatchModelGenerator2D : public ModelGenerator
   static Go::SplineSurface getSubPatch(const Go::SplineSurface* srf,
       const size_t startu, const size_t numcoefsu, const int orderu,
       const size_t startv, const size_t numcoefsv, const int orderv);
+  //! \brief Establishes mixed subdivision bases.
+  //! \param sim Simulator object holding MixedType information
+  bool establishSubdivisionBases (SIMinput& sim) const;
 
 protected:
   //! \brief Generates the G2 description of the geometry.
@@ -145,6 +165,9 @@ class MultiPatchModelGenerator3D : public ModelGenerator
       const size_t startu, const size_t numcoefsu, const int orderu,
       const size_t startv, const size_t numcoefsv, const int orderv,
       const size_t startw, const size_t numcoefsw, const int orderw);
+  //! \brief Establishes mixed subdivision bases.
+  //! \param sim Simulator object holding MixedType information
+  bool establishSubdivisionBases (SIMinput& sim) const;
 
 protected:
   //! \brief Generates the G2 description of the geometry.

Apps/Common/Test/TestMultiPatchModelGenerator.C

@@ -13,6 +13,8 @@
 #include "ASMs1D.h"
 #include "ASMs2D.h"
 #include "ASMs3D.h"
+#include "ASMmxBase.h"
+#include "SAM.h"
 #include "Functions.h"
 #include "IFEM.h"
 #include "MultiPatchModelGenerator.h"
@@ -200,6 +202,45 @@ TEST(TestMultiPatchModelGenerator2D, InnerPatches)
 }
 
 
+TEST(TestMultiPatchModelGenerator2D, SubdivisionsMixed)
+{
+  ASMmxBase::Type = ASMmxBase::FULL_CONT_RAISE_BASIS1;
+  SIMMultiPatchModelGen<SIM2D> sim({1,1});
+  ASSERT_TRUE(sim.read("refdata/modelgen2d_subdivision_mixed.xinp"));
+  const SIM2D::PatchVec& model = sim.getFEModel();
+  sim.preprocess();
+
+  std::vector<std::vector<int>> mlgn =
+        {{1,2,3,4,5,6,7,8,9,10,11,12,13},
+         {2,3,14,5,6,15,8,9,16,11,17,13,18}};
+
+  std::vector<std::vector<int>> mlgeq1 =
+        {{0,1,2,0,3,4,0,5,6,7,8,9,10},
+         {1,2,11,3,4,12,5,6,13,8,14,10,15}};
+
+  for (int i=0; i<2; i++) {
+    check_vector_int_equals(mlgn[i], model[i]->getMyNodeNums());
+    for (size_t n = 0; n < mlgeq1[i].size(); ++n)
+      ASSERT_EQ(sim.getSAM()->getEquation(mlgn[i][n], 1), mlgeq1[i][n]);
+  }
+
+  SIMMultiPatchModelGen<SIM2D> sim2({2,1});
+  ASSERT_TRUE(sim2.read("refdata/modelgen2d_subdivision_mixed.xinp"));
+  sim2.preprocess();
+
+  std::vector<std::vector<int>> mlgeq2 =
+        {{0,1,2,3,4,5,0,6,7,8,9,10,0,11,12,13,14,15,16,17,18,19},
+         {2,3,4,5,20,21,7,8,9,10,22,23,12,13,14,15,24,25,17,26,19,27}};
+
+  for (int i = 0; i < 2; ++i) {
+    size_t k = 0;
+    for (size_t n = 0; n < mlgn[i].size(); ++n)
+      for(size_t d = 1; d <= (n < mlgn[i].size()-4 ? 2 : 1); ++d)
+        ASSERT_EQ(sim2.getSAM()->getEquation(mlgn[i][n], d), mlgeq2[i][k++]);
+  }
+}
+
+
 TEST_P(TestMultiPatchModelGenerator3D, Generate)
 {
   TiXmlDocument doc;
@@ -270,6 +311,46 @@ TEST(TestMultiPatchModelGenerator3D, Subdivisions)
 }
 
 
+TEST(TestMultiPatchModelGenerator3D, SubdivisionsMixed)
+{
+  ASMmxBase::Type = ASMmxBase::FULL_CONT_RAISE_BASIS1;
+  SIMMultiPatchModelGen<SIM3D> sim({1,1});
+  ASSERT_TRUE(sim.read("refdata/modelgen3d_subdivision_mixed.xinp"));
+  const SIM3D::PatchVec& model = sim.getFEModel();
+  sim.preprocess();
+
+  std::vector<std::vector<int>> mlgn =
+        {{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35},
+         {2,3,36,5,6,37,8,9,38,11,12,39,14,15,40,17,18,41,20,21,42,23,24,43,26,27,44,29,45,31,46,33,47,35,48},
+         {4,5,6,7,8,9,49,50,51,13,14,15,16,17,18,52,53,54,22,23,24,25,26,27,55,56,57,30,31,58,59,34,35,60,61},
+         {5,6,37,8,9,38,50,51,62,14,15,40,17,18,41,53,54,63,23,24,43,26,27,44,56,57,64,31,46,59,65,35,48,61,66},
+         {10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,67,68,69,70,71,72,73,74,75,32,33,34,35,76,77,78,79},
+         {11,12,39,14,15,40,17,18,41,20,21,42,23,24,43,26,27,44,68,69,80,71,72,81,74,75,82,33,47,35,48,77,83,79,84},
+         {13,14,15,16,17,18,52,53,54,22,23,24,25,26,27,55,56,57,70,71,72,73,74,75,85,86,87,34,35,60,61,78,79,88,89},
+         {14,15,40,17,18,41,53,54,63,23,24,43,26,27,44,56,57,64,71,72,81,74,75,82,86,87,90,35,48,61,66,79,84,89,91}};
+
+  for (int i=0; i<8; i++)
+    check_vector_int_equals(mlgn[i], model[i]->getMyNodeNums());
+}
+
+
+TEST(TestMultiPatchModelGenerator2D, SubdivisionsPeriodic)
+{
+  ASMmxBase::Type = ASMmxBase::FULL_CONT_RAISE_BASIS1;
+  SIMMultiPatchModelGen<SIM2D> sim({1,1});
+  ASSERT_TRUE(sim.read("refdata/modelgen2d_subdivision_periodic.xinp"));
+  const SIM3D::PatchVec& model = sim.getFEModel();
+  sim.preprocess();
+
+  std::vector<std::vector<int>> mlgn =
+        {{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23},
+         {4,5,1,2,9,10,6,7,14,15,11,12,19,24,16,23,25,20}};
+
+  for (int i=0; i<2; i++)
+    check_vector_int_equals(mlgn[i], model[i]->getMyNodeNums());
+}
+
+
 struct SubPatchTest {
   int n; // number of coefs in total model
   int p; // polynomial order of basis
@@ -592,6 +673,26 @@ INSTANTIATE_TEST_CASE_P(TestGetSubPatch2D,
                         testing::ValuesIn(SubPatch2D));
 
 
+TEST(TestMultiPatchModelGenerator, ExtendedBasis2D)
+{
+  std::vector<double> knots {0.0, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0};
+  Go::BsplineBasis basis(2, knots.begin(), knots.end());
+  Go::BsplineBasis nBasis = extendedBasis(basis);
+  std::vector<double> nknots {0.0, 0.0, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0};
+  check_vector_double_near(nBasis.getKnots(),nknots);
+
+  Go::BsplineBasis basis2(2, knots.begin(), knots.begin()+4);
+  nBasis = extendedBasis(basis2);
+  nknots = {0.0, 0.0, 0.0, 0.2, 0.4, 0.6};
+  check_vector_double_near(nBasis.getKnots(),nknots);
+
+  Go::BsplineBasis basis3(2, knots.begin()+4, knots.end());
+  nBasis = extendedBasis(basis3);
+  nknots = {0.4, 0.6, 0.8, 1.0, 1.0, 1.0};
+  check_vector_double_near(nBasis.getKnots(),nknots);
+}
+
+
 class TestGetSubPatch3D :
   public testing::Test,
   public testing::WithParamInterface<SubPatchTest>

Apps/Common/Test/refdata/modelgen2d_subdivision_mixed.xinp

@@ -0,0 +1,11 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<simulation>
+  <geometry dim="2" sets="true">
+    <subdivision nx="2" ny="1">
+      <refine type="uniform" patch="1" u="1" v="0"/>
+    </subdivision>
+  </geometry>
+  <boundaryconditions>
+    <dirichlet set="Edge1" basis="1" comp="1"/>
+  </boundaryconditions>
+</simulation>

Apps/Common/Test/refdata/modelgen2d_subdivision_periodic.xinp

@@ -0,0 +1,11 @@
+<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
+<simulation>
+  <geometry dim="2" periodic_x="12" sets="true">
+    <subdivision nx="2">
+      <refine patch="1" u="4" v="0"/>
+    </subdivision>
+  </geometry>
+  <boundaryconditions>
+    <dirichlet set="Edge1" basis="1" comp="1"/>
+  </boundaryconditions>
+</simulation>

Apps/Common/Test/refdata/modelgen3d_subdivision_mixed.xinp

@@ -0,0 +1,12 @@
+<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
+<simulation>
+  <geometry sets="true">
+    <partitioning procs="4" nperproc="2"/>
+    <subdivision nx="2" ny="2" nz="2">
+      <refine patch="1" u="1" v="1" w="1"/>
+    </subdivision>
+  </geometry>
+  <boundaryconditions>
+    <dirichlet set="Edge1" basis="1" comp="1"/>
+  </boundaryconditions>
+</simulation>

src/ASM/ASMs2D.C

@@ -632,6 +632,19 @@ bool ASMs2D::connectBasis (int edge, ASMs2D& neighbor, int nedge, bool revers,
   for (int& it : masterNodes)
     it += master;
 
+  int n1, n2;
+  if (!neighbor.getSize(n1,n2,basis)) return false;
+  std::cout << "\tmaster coords:";
+  for (int& i : masterNodes)
+    std::cout << " " << neighbor.getCoord(i)[nedge<3 ? 0 : 1];
+  std::cout << std::endl;
+
+  if (!this->getSize(n1,n2,basis)) return false;
+  std::cout << "\tslave coords:";
+  for (int& i : slaveNodes)
+    std::cout << " " << this->getCoord(i)[nedge<3 ? 0 : 1];
+  std::cout << std::endl;
+
   if (masterNodes.size() != slaveNodes.size())
   {
     std::cerr <<" *** ASMs2D::connectBasis: Non-matching edges, sizes "
@@ -708,15 +721,28 @@ void ASMs2D::constrainEdge (int dir, bool open, int dof, int code, char basis)
   else if (code < 0)
     bcode = -code;
 
+  int i1 = 2;
+  int i2 = 2;
+  int i1e = n1;
+  int i2e = n2;
+  std::vector<int> mult;
   switch (dir)
     {
     case  1: // Right edge (positive I-direction)
       node += n1-1;
     case -1: // Left edge (negative I-direction)
-      if (!open)
-	this->prescribe(node,dof,bcode);
-      node += n1;
-      for (int i2 = 2; i2 < n2; i2++, node += n1)
+      this->getBasis(basis)->basis_v().knotMultiplicities(mult);
+      if (!open) {
+        if (mult[0] == this->getBasis(basis)->basis_v().order()) {
+          this->prescribe(node,dof,bcode);
+          node += n1;
+        } else
+          i2 = 1;
+        if (mult.back() != this->getBasis(basis)->basis_v().order())
+          i2e = n2+1;
+      } else
+        node += n1;
+      for ( ; i2 < i2e; i2++, node += n1)
       {
 	// If the Dirichlet condition is to be projected, add this node to
 	// the set of nodes to receive prescribed value from the projection
@@ -725,16 +751,25 @@ void ASMs2D::constrainEdge (int dir, bool open, int dof, int code, char basis)
 	  dirich.back().nodes.push_back(std::make_pair(i2,node));
       }
       if (!open)
-	this->prescribe(node,dof,bcode);
+        if (mult.back() == this->getBasis(basis)->basis_v().order())
+          this->prescribe(node,dof,bcode);
       break;
 
     case  2: // Back edge (positive J-direction)
       node += n1*(n2-1);
     case -2: // Front edge (negative J-direction)
-      if (!open)
-	this->prescribe(node,dof,bcode);
-      node++;
-      for (int i1 = 2; i1 < n1; i1++, node++)
+      this->getBasis(basis)->basis_u().knotMultiplicities(mult);
+      if (!open) {
+        if (mult[0] == this->getBasis(basis)->basis_u().order()) {
+          this->prescribe(node,dof,bcode);
+          node++;
+        } else
+          i1 = 1;
+        if (mult.back() != this->getBasis(basis)->basis_u().order())
+          i1e = n1+1;
+      } else
+        node++;
+      for ( ; i1 < i1e; i1++, node++)
       {
 	// If the Dirichlet condition is to be projected, add this node to
 	// the set of nodes to receive prescribed value from the projection
@@ -743,7 +778,8 @@ void ASMs2D::constrainEdge (int dir, bool open, int dof, int code, char basis)
 	  dirich.back().nodes.push_back(std::make_pair(i1,node));
       }
       if (!open)
-	this->prescribe(node,dof,bcode);
+        if (mult.back() == this->getBasis(basis)->basis_u().order())
+          this->prescribe(node,dof,bcode);
       break;
     }
 

src/ASM/ASMs2D.h

@@ -115,6 +115,8 @@ class ASMs2D : public ASMstruct, public ASM2D
   //! \brief Returns the spline curve representing a boundary of this patch.
   //! \param[in] dir Parameter direction defining which boundary to return
   virtual Go::SplineCurve* getBoundary(int dir, int = 1);
+  //! \brief Returns the number of bases.
+  virtual size_t getNoBasis() const { return 1; }
   //! \brief Returns the spline surface representing the basis of this patch.
   virtual Go::SplineSurface* getBasis(int = 1) const { return surf; }
   //! \brief Copies the parameter domain from the \a other patch.

src/ASM/ASMs2Dmx.h

@@ -218,6 +218,10 @@ class ASMs2Dmx : public ASMs2D, private ASMmxBase
   virtual void getBoundaryNodes(int lIndex, IntVec& nodes, int basis,
                                 int thick = 1, bool local = false) const;
 
+  //! \brief Manually set the bases for special purposes (e.g. subdivision).
+  //! \param[in] inp The custom constructed basis
+  void setBases(std::vector<std::shared_ptr<Go::SplineSurface>>& inp) { m_basis = inp; }
+
 protected:
   std::vector<std::shared_ptr<Go::SplineSurface>> m_basis; //!< Vector of bases
 };

src/ASM/ASMs3Dmx.h

@@ -206,6 +206,10 @@ class ASMs3Dmx : public ASMs3D, private ASMmxBase
   //! \param[out] n3 Number of nodes in third (w) direction
   //! \param[in] basis Which basis to return size parameters for
   virtual bool getSize(int& n1, int& n2, int& n3, int basis = 0) const;
+
+  //! \brief Manually set the bases for special purposes (e.g. subdivision).
+  //! \param[in] inp The custom constructed basis
+  void setBases(std::vector<std::shared_ptr<Go::SplineVolume>>& inp) { m_basis = inp; }
 protected:
   //! \brief Returns the volume in the parameter space for an element.
   //! \param[in] iel Element index

src/SIM/SIM1D.C

@@ -63,7 +63,8 @@ bool SIM1D::addConnection (int master, int slave, int mIdx, int sIdx,
     if (dim != 0) return false;
 
     IFEM::cout <<"\tConnecting P"<< slave <<" V"<< sIdx
-               <<" to P"<< master <<" V"<< mIdx << std::endl;
+               <<" to P"<< master <<" V"<< mIdx
+               <<" basis="<< basis <<" thick="<< thick << std::endl;
 
     ASMs1D* spch = static_cast<ASMs1D*>(myModel[lslave-1]);
     ASMs1D* mpch = static_cast<ASMs1D*>(myModel[lmaster-1]);

src/SIM/SIM2D.C

@@ -86,7 +86,8 @@ bool SIM2D::addConnection (int master, int slave, int mIdx,
 
     IFEM::cout <<"\tConnecting P"<< slave <<" E"<< sIdx
                <<" to P"<< master <<" E"<< mIdx
-               <<" reversed? "<< orient << std::endl;
+               <<" orient="<< orient <<" basis="<< basis
+               <<" thick="<< thick << std::endl;
 
     ASMs2D* spch = static_cast<ASMs2D*>(myModel[lslave-1]);
     ASMs2D* mpch = static_cast<ASMs2D*>(myModel[lmaster-1]);

src/SIM/SIM3D.C

@@ -64,7 +64,8 @@ bool SIM3D::addConnection (int master, int slave, int mIdx,
 
     IFEM::cout <<"\tConnecting P"<< slave <<" F"<< sIdx
                <<" to P"<< master <<" F"<< mIdx
-               <<" orient "<< orient << std::endl;
+               <<" orient="<< orient <<" basis="<< basis
+               <<" thick="<< thick << std::endl;
 
     ASMs3D* spch = static_cast<ASMs3D*>(myModel[lslave-1]);
     ASMs3D* mpch = static_cast<ASMs3D*>(myModel[lmaster-1]);