Migrate from dune-fufem to dune-functions bases

Migrate from dune-fufem to dune-functions bases.

test/geodesicfeassemblerwrappertest.cc

@@ -65,7 +65,7 @@ typedef std::vector<Rotation<double,dim>> RotationVector;
 typedef MultiTypeBlockVector<DisplacementVector, RotationVector> Vector;
 typedef MultiTypeBlockVector<std::vector<FieldVector<ValueType,dim> >, std::vector<Rotation<ValueType,dim>>> DiffVector; //Equivalent Type used for differentiations
 const int dimCR = Rotation<double,dim>::TangentVector::dimension; //dimCorrectionRotation = Dimension of the correction for rotations
-typedef MultiTypeBlockVector<std::vector<FieldVector<double,dim> >, std::vector<FieldVector<double,dimCR>>> CorrectionType;
+typedef MultiTypeBlockVector<BlockVector<FieldVector<double,dim> >, BlockVector<FieldVector<double,dimCR>>> CorrectionType;
 
 typedef MultiTypeBlockVector<Dune::BCRSMatrix<Dune::FieldMatrix<double,dim,dim>>,  Dune::BCRSMatrix<Dune::FieldMatrix<double,dim,dimCR>>> MatrixRow0;
 typedef MultiTypeBlockVector<Dune::BCRSMatrix<Dune::FieldMatrix<double,dimCR,dim>>, Dune::BCRSMatrix<Dune::FieldMatrix<double,dimCR,dimCR>>> MatrixRow1;
@@ -74,8 +74,8 @@ typedef MultiTypeBlockMatrix<MatrixRow0,MatrixRow1> MatrixType;
 //Types for the Non-mixed space
 typedef RigidBodyMotion<double, dim> RBM;
 const static int blocksize = RBM::TangentVector::dimension;
-typedef Dune::BlockVector<Dune::FieldVector<double, blocksize> > CorrectionTypeNonMixed;
-typedef Dune::BCRSMatrix<Dune::FieldMatrix<double, blocksize, blocksize> > MatrixTypeNonMixed;
+typedef Dune::BlockVector<Dune::FieldVector<double, blocksize> > CorrectionTypeWrapped;
+typedef Dune::BCRSMatrix<Dune::FieldMatrix<double, blocksize, blocksize> > MatrixTypeWrapped;
   
 
 int main (int argc, char *argv[])
@@ -84,20 +84,21 @@ int main (int argc, char *argv[])
 
   FieldVector<double,dim> lower(0), upper(1);
 
-  // ///////////////////////////////////////
-  //    Create the grid
-  // ///////////////////////////////////////
-  using GridType = UGGrid<dim>;
-  auto grid = Dune::StructuredGridFactory<GridType>::createCubeGrid({0,0}, {1,1}, {3,3});
   std::function<bool(Dune::FieldVector<double,dim>)> isSurfaceShell = [](Dune::FieldVector<double,dim> x) {
     return x[2] > 0.99;
   };
-  int refine = 3;
+
+  /////////////////////////////////////////////////////////////////////////
+  //    Create the grid
+  /////////////////////////////////////////////////////////////////////////
+  using GridType = UGGrid<dim>;
+  auto grid = StructuredGridFactory<GridType>::createCubeGrid({0,0,0}, {1,1,1}, {2,1,2});
+  int refine = 1;
   while(refine > 0) {
     for (auto&& e : elements(grid->leafGridView())){
       bool refine = false;
       for (int i = 0; i < e.geometry().corners(); i++) {
-          refine = refine || isSurfaceShell(e.geometry().corner(i));
+        refine = refine || isSurfaceShell(e.geometry().corner(i));
       }
       grid->mark(refine ? 1 : 0, e);
     }
@@ -110,7 +111,7 @@ int main (int argc, char *argv[])
   GridView gridView = grid->leafGridView();
 
   /////////////////////////////////////////////////////////////////////////
-  //  Create a composite basis and the MixedGFEAssembler with its wrapper
+  //  Create a composite basis
   /////////////////////////////////////////////////////////////////////////
 
   using namespace Dune::Functions::BasisFactory;
@@ -126,20 +127,13 @@ int main (int argc, char *argv[])
       )
   ));
 
-  auto deformationPowerBasis = makeBasis(
-    gridView,
-    power<dim>(
-        lagrange<displacementOrder>()
-  ));
   using CompositeBasis = decltype(compositeBasis);
   using LocalView = typename CompositeBasis::LocalView;
 
+  /////////////////////////////////////////////////////////////////////////
+  //  Create the energy functions with their parameters
+  /////////////////////////////////////////////////////////////////////////
 
-  typedef Dune::Functions::LagrangeBasis<GridView,displacementOrder> DeformationFEBasis;
-  typedef Dune::Functions::LagrangeBasis<GridView,rotationOrder> OrientationFEBasis;
-  DeformationFEBasis deformationFEBasis(gridView);
-  OrientationFEBasis orientationFEBasis(gridView);
-  
   ParameterTree parameters;
   //St-Venant-Kirchhoff-Parameters
   parameters["mu"] = "2.7191e+4";
@@ -148,13 +142,12 @@ int main (int argc, char *argv[])
   auto elasticDensity = std::make_shared<Elasticity::StVenantKirchhoffDensity<dim,ValueType>>(parameters);
   auto elasticEnergy = std::make_shared<GFE::LocalIntegralEnergy<LocalView, ValueType>>(elasticDensity);
 
-  //prepare stuff for the surfacecosseratenergy...
   //Surface-Cosserat-Energy-Parameters
-  parameters["mu_c_"] = "0";
-  parameters["L_c_"] = "0.001";
-  parameters["b1_"] = "1";
-  parameters["b2_"] = "1";
-  parameters["b3_"] = "1";
+  parameters["mu_c"] = "0";
+  parameters["L_c"] = "0.001";
+  parameters["b1"] = "1";
+  parameters["b2"] = "1";
+  parameters["b3"] = "1";
 
   std::vector<UnitVector<double,dim> > vertexNormals(gridView.size(dim));
   Dune::FieldVector<double,dim> vertexNormalRaw = {0,0,1};
@@ -209,13 +202,19 @@ int main (int argc, char *argv[])
                     Rotation<double,dim> > mixedAssembler(compositeBasis, &mixedLocalGFEADOLCStiffness);
 
   typedef RigidBodyMotion<double, dim> RBM;
+  typedef Dune::Functions::LagrangeBasis<GridView,displacementOrder> DeformationFEBasis;
+  DeformationFEBasis deformationFEBasis(gridView);
   typedef GeodesicFEAssemblerWrapper<CompositeBasis, DeformationFEBasis, RBM, RealTuple<double, dim>, Rotation<double,dim>> GFEAssemblerWrapper;
   GFEAssemblerWrapper assembler(&mixedAssembler, deformationFEBasis);
-
-  //////////////////////////////////////////////////////////////
+  
+  /////////////////////////////////////////////////////////////////////////
   //  Prepare the iterate x where we want to assemble
-  //////////////////////////////////////////////////////////////  
-
+  /////////////////////////////////////////////////////////////////////////  
+  auto deformationPowerBasis = makeBasis(
+    gridView,
+    power<dim>(
+        lagrange<displacementOrder>()
+  ));
   BlockVector<FieldVector<double,dim> > identity(compositeBasis.size({0}));
   Dune::Functions::interpolate(deformationPowerBasis, identity, [](FieldVector<double,dim> x){ return x; });
 
@@ -233,16 +232,63 @@ int main (int argc, char *argv[])
   }
 
   //////////////////////////////////////////////////////////////
-  //  Compute the energy, assemble and compare!
+  //  Compute the energy, assemble the  and compare!
   //////////////////////////////////////////////////////////////
   
-  CorrectionTypeNonMixed gradient;
-  MatrixTypeNonMixed hessianMatrix;
+  CorrectionTypeWrapped gradient;
+  MatrixTypeWrapped hessianMatrix;
   double energy = assembler.computeEnergy(xRBM);
   assembler.assembleGradientAndHessian(xRBM, gradient, hessianMatrix, true);
+  double gradientTwoNorm = gradient.two_norm();
+  double gradientInfinityNorm = gradient.infinity_norm();
+  double matrixFrobeniusNorm = hessianMatrix.frobenius_norm();
 
   CorrectionType gradientMixed;
+  gradientMixed[_0].resize(x[_0].size());
+  gradientMixed[_1].resize(x[_1].size());
   MatrixType hessianMatrixMixed;
-  //double energyMixed = mixedAssembler.computeEnergy(x);
-  //mixedAssembler.assembleGradientAndHessian(x, gradientMixed, hessianMatrixMixed, true);
+  double energyMixed = mixedAssembler.computeEnergy(x[_0], x[_1]);
+  mixedAssembler.assembleGradientAndHessian(x[_0], x[_1], gradientMixed[_0], gradientMixed[_1], hessianMatrixMixed, true);
+  double gradientMixedTwoNorm = gradientMixed.two_norm();
+  double gradientMixedInfinityNorm = gradientMixed.infinity_norm();
+  double matrixMixedFrobeniusNorm = hessianMatrixMixed.frobenius_norm();
+
+  //These values were taken from assembling using a fufem-basis
+  double expectedEnergy = 0.103846154;
+  double expectedGradientTwoNorm = 0.526624705;
+  double expectedGradientInfinityNorm = 0.230769231;
+  double expectedMatrixFrobeniusNorm = 1718346.31;
+
+  if ( std::abs(energy - expectedEnergy)/expectedEnergy > 1e-4 || std::abs(energy - energyMixed)/energyMixed > 1e-4)
+  {
+      std::cerr << std::setprecision(9);
+      std::cerr << "The energy calculated by the GeodesicFEAssemblerWrapper is " << energy << " but "
+                << expectedEnergy << " was expected!" << std::endl;
+      std::cerr << "The energy calculated by the MixedGFEAssembler is " << energyMixed << "!" << std::endl;
+      return 1;
+  }
+  if ( std::abs(gradientTwoNorm - expectedGradientTwoNorm)/expectedGradientTwoNorm > 1e-4 ||
+       std::abs(gradientInfinityNorm - expectedGradientInfinityNorm)/expectedGradientInfinityNorm > 1e-8 ||
+       std::abs(gradientTwoNorm - gradientMixedTwoNorm)/gradientMixedTwoNorm > 1e-4 ||
+       std::abs(gradientInfinityNorm - gradientMixedInfinityNorm)/gradientMixedInfinityNorm > 1e-8)
+  {
+      std::cerr << std::setprecision(9);
+      std::cerr << "The gradient infinity norm calculated by the GeodesicFEAssemblerWrapper is " << gradientInfinityNorm << " but "
+                << expectedGradientInfinityNorm << " was expected!" << std::endl;
+      std::cerr << "The gradient infinity norm calculated by the MixedGFEAssembler is " << gradientMixedInfinityNorm << "!" << std::endl;
+      std::cerr << "The gradient norm calculated by the GeodesicFEAssemblerWrapper is " << gradientTwoNorm << " but "
+                << expectedGradientTwoNorm << " was expected!" << std::endl;
+      std::cerr << "The gradient norm calculated by the MixedGFEAssembler is " << gradientMixedTwoNorm << "!" << std::endl;
+      return 1;
+  }
+
+  if ( std::abs(matrixFrobeniusNorm - expectedMatrixFrobeniusNorm)/expectedMatrixFrobeniusNorm > 1e-4 ||
+       std::abs(matrixFrobeniusNorm - matrixMixedFrobeniusNorm)/matrixMixedFrobeniusNorm > 1e-4 )
+  {
+      std::cerr << std::setprecision(9);
+      std::cerr << "The matrix norm calculated by the GeodesicFEAssemblerWrapper is " << matrixFrobeniusNorm << " but "
+                << expectedMatrixFrobeniusNorm << " was expected!" << std::endl;
+      std::cerr << "The matrix norm calculated by the MixedGFEAssembler is " << matrixMixedFrobeniusNorm << "!" << std::endl;
+      return 1;
+  }
 }