changed: move 3D mixed establishment code to ASMmxBase

mirrors the 2D code

src/ASM/ASMmxBase.C

@@ -14,6 +14,8 @@
 #include "ASMmxBase.h"
 #include "GoTools/geometry/SplineSurface.h"
 #include "GoTools/geometry/SurfaceInterpolator.h"
+#include "GoTools/trivariate/SplineVolume.h"
+#include "GoTools/trivariate/VolumeInterpolator.h"
 
 char ASMmxBase::geoBasis            = 2;
 ASMmxBase::MixedType ASMmxBase::Type = ASMmxBase::FULL_CONT_RAISE_BASIS1;
@@ -184,3 +186,59 @@ ASMmxBase::SurfaceVec ASMmxBase::establishBases(Go::SplineSurface* surf,
 
   return result;
 }
+
+
+ASMmxBase::VolumeVec ASMmxBase::establishBases(Go::SplineVolume* svol,
+                                               MixedType type)
+{
+  VolumeVec result(2);
+  // With mixed methods we need two separate spline spaces
+  if (type == FULL_CONT_RAISE_BASIS1 || type == FULL_CONT_RAISE_BASIS2)
+  {
+    // basis1 should be one degree higher than basis2 and C^p-1 continuous
+    int ndim = svol->dimension();
+    Go::BsplineBasis b1 = svol->basis(0).extendedBasis(svol->order(0)+1);
+    Go::BsplineBasis b2 = svol->basis(1).extendedBasis(svol->order(1)+1);
+    Go::BsplineBasis b3 = svol->basis(2).extendedBasis(svol->order(2)+1);
+    /* To lower order and regularity this can be used instead
+    std::vector<double>::const_iterator first = ++surf->basis(0).begin();
+    std::vector<double>::const_iterator last  = --surf->basis(0).end();
+    Go::BsplineBasis b1 = Go::BsplineBasis(surf->order_u()-1,first,last);
+    first =  ++surf->basis(1).begin();
+    last  =  --surf->basis(1).end();
+    Go::BsplineBasis b2 = Go::BsplineBasis(surf->order_v()-1,first,last);
+    */
+
+    // Compute parameter values of the Greville points
+    size_t i;
+    RealArray ug(b1.numCoefs()), vg(b2.numCoefs()), wg(b3.numCoefs());
+    for (i = 0; i < ug.size(); i++)
+      ug[i] = b1.grevilleParameter(i);
+    for (i = 0; i < vg.size(); i++)
+      vg[i] = b2.grevilleParameter(i);
+    for (i = 0; i < wg.size(); i++)
+      wg[i] = b3.grevilleParameter(i);
+
+    RealArray XYZ(ndim*ug.size()*vg.size()*wg.size());
+    // Evaluate the spline surface at all points
+    svol->gridEvaluator(XYZ,ug,vg,wg);
+
+    // Project the coordinates onto the new basis (the 2nd XYZ is dummy here)
+    result[0].reset(Go::VolumeInterpolator::regularInterpolation(b1,b2,b3,
+                                                                 ug,vg,wg,XYZ,ndim,
+                                                                 false,XYZ));
+  }
+  else if (type == REDUCED_CONT_RAISE_BASIS1 || type == REDUCED_CONT_RAISE_BASIS2)
+  {
+    // Order-elevate basis1 such that it is of one degree higher than basis2
+    // but only C^p-2 continuous
+    result[0].reset(new Go::SplineVolume(*svol));
+    result[0]->raiseOrder(1,1,1);
+  }
+  result[1].reset(new Go::SplineVolume(*svol));
+
+  if (type == FULL_CONT_RAISE_BASIS2 || type == REDUCED_CONT_RAISE_BASIS2)
+    std::swap(result[0], result[1]);
+
+  return result;
+}

src/ASM/ASMmxBase.h

@@ -19,6 +19,7 @@
 
 namespace Go {
   class SplineSurface;
+  class SplineVolume;
 }
 
 
@@ -74,6 +75,7 @@ public:
   static MixedType Type; //!< Type of mixed formulation used
 
   typedef std::vector<std::shared_ptr<Go::SplineSurface>> SurfaceVec; //!< Convenience type
+  typedef std::vector<std::shared_ptr<Go::SplineVolume>> VolumeVec; //!< Convenience type
 
   //! \brief Establish mixed bases
   //! \param[in] surf The base basis to use.
@@ -81,6 +83,12 @@ public:
   //! \return Vector with bases.
   static SurfaceVec establishBases(Go::SplineSurface* surf, MixedType type);
 
+  //! \brief Establish mixed bases
+  //! \param[in] vol The base basis to use.
+  //! \param[in] type The type of bases to establish.
+  //! \return Vector with bases.
+  static VolumeVec establishBases(Go::SplineVolume* vol, MixedType type);
+
 private:
   std::vector<int> MADOF; //!< Matrix of accumulated DOFs for this patch
 

src/ASM/ASMs3Dmx.C

@@ -175,73 +175,7 @@ bool ASMs3Dmx::generateFEMTopology ()
   if (!svol) return false;
 
   if (m_basis.empty())
-  {
+    m_basis = ASMmxBase::establishBases(svol, ASMmxBase::Type);
-    m_basis.resize(2);
-    // With mixed methods we need two separate spline spaces
-    if (ASMmxBase::Type == FULL_CONT_RAISE_BASIS1 ||
-        ASMmxBase::Type == FULL_CONT_RAISE_BASIS2)
-    {
-      // basis1 should be one degree higher than basis2 and C^p-1 continuous
-      int ndim = svol->dimension();
-      Go::BsplineBasis b1 = svol->basis(0).extendedBasis(svol->order(0)+1);
-      Go::BsplineBasis b2 = svol->basis(1).extendedBasis(svol->order(1)+1);
-      Go::BsplineBasis b3 = svol->basis(2).extendedBasis(svol->order(2)+1);
-      /* To lower order and regularity this can be used instead
-      std::vector<double>::const_iterator first =  ++svol->basis(0).begin();
-      std::vector<double>::const_iterator last  =  --svol->basis(0).end();
-      Go::BsplineBasis b1 = Go::BsplineBasis(svol->order(0)-1,first,last);
-      first =  ++svol->basis(1).begin();
-      last  =  --svol->basis(1).end();
-      Go::BsplineBasis b2 = Go::BsplineBasis(svol->order(1)-1,first,last);
-      first =  ++svol->basis(2).begin();
-      last  =  --svol->basis(2).end();
-      Go::BsplineBasis b3 = Go::BsplineBasis(svol->order(2)-1,first,last);
-      */
-
-      // Note: Currently this is implemented for non-rational splines only.
-      // TODO: Ask the splines people how to fix this properly, that is, how
-      // may we obtain the correct weights for basis1 when *svol is a NURBS?
-      if (svol->rational())
-	std::cout <<"WARNING: The geometry basis is rational (using NURBS).\n"
-		  <<"         The basis for the unknown fields of one degree"
-		  <<" higher will however be non-rational.\n"
-		  <<"         This may affect accuracy.\n"<< std::endl;
-
-      // Compute parameter values of the Greville points of the new basis
-      size_t i;
-      RealArray ug(b1.numCoefs()), vg(b2.numCoefs()), wg(b3.numCoefs());
-      for (i = 0; i < ug.size(); i++)
-	ug[i] = b1.grevilleParameter(i);
-      for (i = 0; i < vg.size(); i++)
-	vg[i] = b2.grevilleParameter(i);
-      for (i = 0; i < wg.size(); i++)
-	wg[i] = b3.grevilleParameter(i);
-
-      // Evaluate the spline volume at all points
-      RealArray XYZ(ndim*ug.size()*vg.size()*wg.size());
-      svol->gridEvaluator(ug,vg,wg,XYZ);
-
-      // Project the coordinates onto the new basis (the 2nd XYZ is dummy here)
-      m_basis[0].reset(Go::VolumeInterpolator::regularInterpolation(b1,b2,b3,
-							            ug,vg,wg,XYZ,ndim,
-                                                                    false,XYZ));
-    }
-    else if (ASMmxBase::Type == REDUCED_CONT_RAISE_BASIS1 ||
-             ASMmxBase::Type == REDUCED_CONT_RAISE_BASIS2)
-    {
-      // Order-elevate basis1 such that it is of one degree higher than basis2
-      // but only C^p-2 continuous
-      m_basis[0].reset(new Go::SplineVolume(*svol));
-      m_basis[0]->raiseOrder(1,1,1);
-    }
-    m_basis[1].reset(new Go::SplineVolume(*svol));
-
-    if (ASMmxBase::Type == FULL_CONT_RAISE_BASIS2 ||
-        ASMmxBase::Type == REDUCED_CONT_RAISE_BASIS2) {
-      std::swap(m_basis[0], m_basis[1]);
-      svol = m_basis[0].get();
-    }
-  }
 
   nb.clear();
   for (auto it : m_basis)