Codipack parallel and mooneyrivlin code

dune/elasticity/assemblers/feassembler.hh

@@ -11,7 +11,7 @@
 
 /** \brief A global FE assembler for problems involving functions that map into non-Euclidean spaces
  */
-template <class Basis, class VectorType>
+template <class Basis, class VectorType, class GradientVectorType>
 class FEAssembler {
 
     typedef typename Basis::GridView GridView;
@@ -32,13 +32,13 @@ protected:
 
     LocalFEStiffness<GridView,
                              typename Basis::LocalFiniteElement,
-                             VectorType>* localStiffness_;
+                             VectorType, GradientVectorType>* localStiffness_;
 
 public:
 
     /** \brief Constructor for a given grid */
     FEAssembler(const Basis& basis,
-                LocalFEStiffness<GridView,typename Basis::LocalFiniteElement, VectorType>* localStiffness)
+                LocalFEStiffness<GridView,typename Basis::LocalFiniteElement, VectorType, GradientVectorType>* localStiffness)
         : basis_(basis),
           localStiffness_(localStiffness)
     {}
@@ -63,15 +63,15 @@ public:
 
 
 
-template <class Basis, class TargetSpace>
-void FEAssembler<Basis,TargetSpace>::
+template <class Basis, class TargetSpace, class GradientVectorType>
+void FEAssembler<Basis,TargetSpace, GradientVectorType>::
 getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
 {
     int n = basis_.size();
 
     nb.resize(n, n);
 
-    for (const auto& element : elements(basis_.getGridView(), Dune::Partitions::interior))
+    for (const auto& element : elements(basis_.getGridView()))
     {
         const typename Basis::LocalFiniteElement& lfe = basis_.getLocalFiniteElement(element);
 
@@ -92,8 +92,8 @@ getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
 
 }
 
-template <class Basis, class VectorType>
-void FEAssembler<Basis,VectorType>::
+template <class Basis, class VectorType, class GradientVectorType>
+void FEAssembler<Basis,VectorType, GradientVectorType>::
 assembleGradientAndHessian(const VectorType& sol,
                            Dune::BlockVector<Dune::FieldVector<double, blocksize> >& gradient,
                            Dune::BCRSMatrix<MatrixBlock>& hessian,
@@ -112,7 +112,7 @@ assembleGradientAndHessian(const VectorType& sol,
     gradient.resize(sol.size());
     gradient = 0;
 
-    for (const auto& element : elements(basis_.getGridView(), Dune::Partitions::interior))
+    for (const auto& element : elements(basis_.getGridView()))
     {
         const int numOfBaseFct = basis_.getLocalFiniteElement(element).localBasis().size();
 
@@ -150,8 +150,8 @@ assembleGradientAndHessian(const VectorType& sol,
 }
 
 
-template <class Basis, class VectorType>
-double FEAssembler<Basis, VectorType>::
+template <class Basis, class VectorType, class GradientVectorType>
+double FEAssembler<Basis, VectorType, GradientVectorType>::
 computeEnergy(const VectorType& sol) const
 {
     double energy = 0;
@@ -174,7 +174,10 @@ computeEnergy(const VectorType& sol) const
 
     }
 
-    return energy;
-
+    auto&& grid = basis_.getGridView().grid();
+    double totalEnergy = grid.comm().sum(energy);
+    std::cout << "rank " << grid.comm().rank() << " gives energy of " << energy << " and total = " << totalEnergy << std::endl;
+    grid.comm().broadcast(&totalEnergy, 1, 0);
+    return totalEnergy;
 }
 #endif