Fixed: Mixed phase-field solver, assign zero dofs for second basis.

Removed unused 1D version, no longer compile due to Dim::nf reference.

Author: kmokstad

CahnHilliard.C

@@ -305,6 +305,13 @@ bool CahnHilliard::evalInt (LocalIntegral& elmInt, const FiniteElement& fe,
 }
 
 
+bool CahnHilliard::evalIntMx (LocalIntegral& elmInt, const MxFiniteElement& fe,
+                              const Vec3& X) const
+{
+  return this->evalInt(elmInt,fe,X);
+}
+
+
 bool CahnHilliard::evalBou (LocalIntegral& elmInt, const FiniteElement& fe,
                             const Vec3& X, const Vec3& normal) const
 {
@@ -363,6 +370,12 @@ NormBase* CahnHilliard::getNormIntegrand (AnaSol* a) const
 }
 
 
+CahnHilliard4::CahnHilliard4 (unsigned short int n) : CahnHilliard(n)
+{
+  scale2nd = 2.0;
+}
+
+
 bool CahnHilliard4::evalInt (LocalIntegral& elmInt, const FiniteElement& fe,
                              const Vec3& X) const
 {
@@ -460,6 +473,14 @@ bool CahnHilliardNorm::evalInt (LocalIntegral& elmInt, const FiniteElement& fe,
 }
 
 
+bool CahnHilliardNorm::evalIntMx (LocalIntegral& elmInt,
+                                  const MxFiniteElement& fe,
+                                  const Vec3& X) const
+{
+  return this->evalInt(elmInt,fe,X);
+}
+
+
 bool CahnHilliardNorm::finalizeElement (LocalIntegral& elmInt)
 {
   if (Lnorm == 0) return true;

CahnHilliard.h

@@ -64,6 +64,14 @@ class CahnHilliard : public IntegrandBase
   virtual bool evalInt(LocalIntegral& elmInt, const FiniteElement& fe,
                        const Vec3& X) const;
 
+  using IntegrandBase::evalIntMx;
+  //! \brief Evaluates the integrand at an interior point.
+  //! \param elmInt The local integral object to receive the contributions
+  //! \param[in] fe Mixed finite element data of current integration point
+  //! \param[in] X Cartesian coordinates of current integration point
+  virtual bool evalIntMx(LocalIntegral& elmInt, const MxFiniteElement& fe,
+                         const Vec3& X) const;
+
   using IntegrandBase::evalSol;
   //! \brief Evaluates the secondary solution at a result point.
   //! \param[out] s Array of solution field values at current point
@@ -156,7 +164,7 @@ class CahnHilliard4 : public CahnHilliard
 {
 public:
   //! \brief The constructor forwards to the parent class constructor.
-  explicit CahnHilliard4(unsigned short int n) : CahnHilliard(n) { scale2nd = 2.0; }
+  explicit CahnHilliard4(unsigned short int n);
   //! \brief Empty destructor.
   virtual ~CahnHilliard4() {}
 
@@ -193,6 +201,14 @@ class CahnHilliardNorm : public NormBase
   virtual bool evalInt(LocalIntegral& elmInt, const FiniteElement& fe,
                        const Vec3& X) const;
 
+  using NormBase::evalIntMx;
+  //! \brief Evaluates the integrand at an interior point.
+  //! \param elmInt The local integral object to receive the contributions
+  //! \param[in] fe Mixed finite element data of current integration point
+  //! \param[in] X Cartesian coordinates of current integration point
+  virtual bool evalIntMx(LocalIntegral& elmInt, const MxFiniteElement& fe,
+                         const Vec3& X) const;
+
   using NormBase::finalizeElement;
   //! \brief Finalizes the element norms after the numerical integration.
   //! \param elmInt The local integral object to receive the contributions

SIMPhaseField.h

@@ -55,7 +55,10 @@ template<class Dim> class SIMPhaseField : public Dim
   {
     Dim::myHeading = "Cahn-Hilliard solver";
     if (gridOwner && gridOwner->createFEMmodel())
+    {
+      Dim::nf.resize(gridOwner->getNoBasis(),0);
       this->clonePatches(gridOwner->getFEModel(),gridOwner->getGlob2LocMap());
+    }
 
     eps_d0 = refTol = 0.0;
     vtfStep = irefine = 0;

main_CahnHilliard.C

@@ -12,7 +12,6 @@
 //==============================================================================
 
 #include "IFEM.h"
-#include "SIM1D.h"
 #include "SIM2D.h"
 #include "SIM3D.h"
 #include "SIMPhaseField.h"
@@ -82,8 +81,6 @@ int main (int argc, char** argv)
       ; // ignore the obsolete option
     else if (!strcmp(argv[i],"-2D"))
       ndim = 2;
-    else if (!strcmp(argv[i],"-1D"))
-      ndim = 1;
     else if (!infile)
       infile = argv[i];
     else
@@ -93,7 +90,7 @@ int main (int argc, char** argv)
   {
     std::cout <<"usage: "<< argv[0]
               <<" <inputfile> [-dense|-spr|-superlu[<nt>]|-samg|-petsc]\n"
-              <<"       [-lag|-spec|-LR] [-1D|-2D] [-nGauss <n>]\n"
+              <<"       [-lag|-spec|-LR] [-2D] [-nGauss <n>]\n"
               <<"       [-vtf <format> [-nviz <nviz>]"
               <<" [-nu <nu>] [-nv <nv>] [-nw <nw>]] [-hdf5]"<< std::endl;
     return 0;
@@ -107,5 +104,5 @@ int main (int argc, char** argv)
   else if (ndim == 2)
     return runSimulator<SIM2D>(infile);
   else
-    return runSimulator<SIM1D>(infile);
+    return 1; // 1D no longer supported
 }