hyperelasticity_functional.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
#include <kernel/base_header.hpp>
#include <kernel/assembly/slip_filter_assembler.hpp>
#include <kernel/assembly/unit_filter_assembler.hpp>
#include <kernel/geometry/export_vtk.hpp>
#include <kernel/geometry/mesh_node.hpp>
#include <kernel/lafem/filter_chain.hpp>
#include <kernel/lafem/filter_sequence.hpp>
#include <kernel/lafem/slip_filter.hpp>
#include <kernel/lafem/unit_filter_blocked.hpp>
#include <kernel/meshopt/mesh_concentration_function.hpp>
#include <kernel/meshopt/mesh_quality_functional.hpp>
#include <kernel/meshopt/rumpf_trafo.hpp>
#include <kernel/util/dist.hpp>
 
#include <map>
 
// This is for the explicit template instantiation below
#ifdef FEAT_EICKT
#include <kernel/meshopt/rumpf_functional.hpp>
#include <kernel/meshopt/rumpf_functionals/p1.hpp>
#include <kernel/meshopt/rumpf_functionals/q1.hpp>
#endif
 
 
namespace FEAT
{
  namespace Meshopt
  {
 
    enum class ScaleComputation
    {
      undefined = 0,
      once_uniform,
      once_cellsize,
      once_concentration,
      current_uniform,
      current_cellsize,
      current_concentration,
      iter_concentration
    };
 
 
    inline std::ostream& operator<<(std::ostream& os, ScaleComputation scale_computation)
    {
      switch(scale_computation)
      {
        case ScaleComputation::undefined:
          return os << "undefined";
        case ScaleComputation::once_uniform:
          return os << "once_uniform";
        case ScaleComputation::once_cellsize:
          return os << "once_cellsize";
        case ScaleComputation::once_concentration:
          return os << "once_concentration";
        case ScaleComputation::current_uniform:
          return os << "current_uniform";
        case ScaleComputation::current_cellsize:
          return os << "current_cellsize";
        case ScaleComputation::current_concentration:
          return os << "current_concentration";
        case ScaleComputation::iter_concentration:
          return os << "iter_concentration";
        default:
          return os << "-unknown-";
      }
    }
 
    inline void operator<<(ScaleComputation& scale_computation, const String& sc_name)
    {
      if(sc_name == "undefined")
        scale_computation = ScaleComputation::undefined;
      else if(sc_name == "once_uniform")
        scale_computation = ScaleComputation::once_uniform;
      else if(sc_name == "once_cellsize")
        scale_computation = ScaleComputation::once_cellsize;
      else if(sc_name == "once_concentration")
        scale_computation = ScaleComputation::once_concentration;
      else if(sc_name == "current_uniform")
        scale_computation = ScaleComputation::current_uniform;
      else if(sc_name == "current_cellsize")
        scale_computation = ScaleComputation::current_cellsize;
      else if(sc_name == "current_concentration")
        scale_computation = ScaleComputation::current_concentration;
      else if(sc_name == "iter_concentration")
        scale_computation = ScaleComputation::iter_concentration;
      else
        XABORTM("Unknown ScaleComputation identifier string " +sc_name);
    }
 
    template
    <
      typename DT_,
      typename IT_,
      typename Trafo_,
      typename CellFunctionalType_,
      typename RefCellTrafo_ = RumpfTrafo<Trafo_, typename Trafo_::MeshType::CoordType>
    >
    class HyperelasticityFunctional :
      public MeshQualityFunctional<typename Trafo_::MeshType>
    {
      public :
        typedef Trafo_ TrafoType;
        typedef typename TrafoType::MeshType MeshType;
        typedef typename MeshType::CoordType CoordType;
        typedef CellFunctionalType_ CellFunctionalType;
 
        typedef RefCellTrafo_ RefCellTrafo;
 
        typedef CoordType DataType;
        typedef Index IndexType;
 
        typedef MeshQualityFunctional<MeshType> BaseClass;
 
        static constexpr int BlockHeight = MeshType::world_dim;
        static constexpr int BlockWidth = MeshType::world_dim;
 
        typedef typename MeshType::ShapeType ShapeType;
 
        typedef LAFEM::DenseVector<CoordType, IndexType> ScalarVectorType;
        typedef LAFEM::DenseVectorBlocked<CoordType, IndexType, MeshType::world_dim> VectorType;
        typedef LAFEM::UnitFilterBlocked<DT_, IT_, MeshType::world_dim> DirichletFilterType;
        typedef LAFEM::FilterSequence<DirichletFilterType> DirichletFilterSequence;
        typedef LAFEM::SlipFilter<DT_, IT_, MeshType::world_dim> SlipFilterType;
        typedef LAFEM::FilterSequence<SlipFilterType> SlipFilterSequence;
        typedef LAFEM::FilterChain<SlipFilterSequence, DirichletFilterSequence> FilterType;
 
        typedef typename Intern::TrafoFE<TrafoType>::Space SpaceType;
 
        typedef LAFEM::DenseVectorBlocked<DT_, IT_, MeshType::world_dim> VectorTypeL;
        typedef LAFEM::DenseVectorBlocked<DT_, IT_, MeshType::world_dim> VectorTypeR;
        typedef VectorTypeR GradientType;
        typedef typename BaseClass::CoordsBufferType CoordsBufferType;
 
        static_assert(std::is_same<ShapeType, typename CellFunctionalType::ShapeType>::value,
        "ShapeTypes of the transformation / cell functional have to agree" );
 
      protected:
        std::shared_ptr<CellFunctionalType> _cell_functional;
        TrafoType& _trafo;
        ScalarVectorType _mu;
        // In the optimal case, every cell in the mesh has the size lambda(cell)
        ScalarVectorType _lambda;
        ScalarVectorType _h;
        std::map<String, std::shared_ptr<Assembly::UnitFilterAssembler<MeshType>>> _dirichlet_asm;
        std::map<String, std::shared_ptr<Assembly::SlipFilterAssembler<TrafoType>>> _slip_asm;
        std::shared_ptr<MeshConcentrationFunctionBase<Trafo_, RefCellTrafo_>> _mesh_conc;
 
      public:
        SpaceType trafo_space;
        std::map<String, DataType*> sync_scalars;
        std::set<VectorTypeR*> sync_vecs;
 
      private:
        const Index _columns;
        const Index _rows;
        ScaleComputation _scale_computation;
 
        DataType _penalty_param;
        DataType _alignment_constraint;
 
        DataType _sum_det;
        DataType _sum_lambda;
        DataType _sum_mu;
 
      public:
        explicit HyperelasticityFunctional(
          Geometry::RootMeshNode<MeshType>* rmn_,
          TrafoType& trafo,
          const std::deque<String>& dirichlet_list,
          const std::deque<String>& slip_list,
          std::shared_ptr<CellFunctionalType_> cell_functional_)
          : BaseClass(rmn_),
          _cell_functional(cell_functional_),
          _trafo(trafo),
          _mu(rmn_->get_mesh()->get_num_entities(ShapeType::dimension), DataType(1)),
          _lambda(rmn_->get_mesh()->get_num_entities(ShapeType::dimension)),
          _h(rmn_->get_mesh()->get_num_entities(ShapeType::dimension)),
          _dirichlet_asm(),
          _slip_asm(),
          _mesh_conc(nullptr),
          trafo_space(trafo),
          sync_scalars(),
          sync_vecs(),
          _columns(trafo_space.get_num_dofs()),
          _rows(trafo_space.get_num_dofs()),
          _scale_computation(ScaleComputation::once_uniform),
          _penalty_param(0),
          _alignment_constraint(0),
          _sum_det(0),
          _sum_lambda(0),
          _sum_mu(0)
          {
 
            XASSERTM(cell_functional_ != nullptr, "cell functional must not be nullptr");
 
            // For every MeshPart specified in the list of Dirichlet boundaries, create a UnitFilterAssembler and
            // insert it into the std::map
            for(auto& it : dirichlet_list)
            {
              // Create empty assembler
              auto new_asm = std::make_shared<Assembly::UnitFilterAssembler<MeshType>>();
 
              // Add the MeshPart to the assembler if it is there. There are legitimate reasons for it NOT to be
              // there, i.e. we are in parallel and our patch is not adjacent to that MeshPart
              auto* mpp = rmn_->find_mesh_part(it);
              if(mpp != nullptr)
              {
                new_asm->add_mesh_part(*mpp);
              }
 
              // Insert into the map
              String identifier(it);
              _dirichlet_asm.emplace(identifier, new_asm);
            }
 
            // For every MeshPart specified in the list of slip boundaries, create a SlipFilterAssembler and
            // insert it into the std::map
            for(auto& it : slip_list)
            {
              // Create empty assembler
              auto new_asm = std::make_shared<Assembly::SlipFilterAssembler<TrafoType>>(this->_trafo);
 
              // Add the MeshPart to the assembler if it is there. There are legitimate reasons for it NOT to be
              // there, i.e. we are in parallel and our patch is not adjacent to that MeshPart
              auto* mpp = rmn_->find_mesh_part(it);
              if(mpp != nullptr)
              {
                new_asm->add_mesh_part(*mpp);
              }
 
              // Insert into the map
              String identifier(it);
              _slip_asm.emplace(identifier, new_asm);
            }
 
            for(Index i(0); i < _mu.size(); ++i)
            {
              _sum_mu += _mu(i);
            }
 
            sync_scalars.emplace("_sum_mu",&_sum_mu);
 
            // Compute desired element size distribution
            _compute_scales_once();
          }
 
        explicit HyperelasticityFunctional(
          Geometry::RootMeshNode<MeshType>* rmn_,
          TrafoType& trafo,
          const std::deque<String>& dirichlet_list,
          const std::deque<String>& slip_list,
          std::shared_ptr<CellFunctionalType_> cell_functional_,
          ScaleComputation scale_computation_,
          std::shared_ptr<MeshConcentrationFunctionBase<Trafo_, RefCellTrafo_>> mesh_conc_ = nullptr,
          DataType penalty_param_ = DataType(0))
          : BaseClass(rmn_),
          _cell_functional(cell_functional_),
          _trafo(trafo),
          _mu(rmn_->get_mesh()->get_num_entities(ShapeType::dimension), DataType(1)),
          _lambda(rmn_->get_mesh()->get_num_entities(ShapeType::dimension)),
          _h(rmn_->get_mesh()->get_num_entities(ShapeType::dimension)),
          _dirichlet_asm(),
          _slip_asm(),
          _mesh_conc(mesh_conc_),
          trafo_space(trafo),
          sync_scalars(),
          sync_vecs(),
          _columns(trafo_space.get_num_dofs()),
          _rows(trafo_space.get_num_dofs()),
          _scale_computation(scale_computation_),
          _penalty_param(penalty_param_),
          _alignment_constraint(0),
          _sum_det(0),
          _sum_lambda(0),
          _sum_mu(0)
          {
 
            XASSERTM(cell_functional_ != nullptr, "cell functional must not be nullptr!\n");
            XASSERTM(_penalty_param >= DataType(0), "penalty_param must be >= 0!\n");
 
            // For every MeshPart specified in the list of Dirichlet boundaries, create a UnitFilterAssembler and
            // insert it into the std::map
            for(auto& it : dirichlet_list)
            {
              // Create empty assembler
              auto new_asm = std::make_shared<Assembly::UnitFilterAssembler<MeshType>>();
 
              // Add the MeshPart to the assembler if it is there. There are legitimate reasons for it NOT to be
              // there, i.e. we are in parallel and our patch is not adjacent to that MeshPart
              auto* mpp = rmn_->find_mesh_part(it);
              if(mpp != nullptr)
              {
                new_asm->add_mesh_part(*mpp);
              }
 
              // Insert into the map
              String identifier(it);
              _dirichlet_asm.emplace(identifier, new_asm);
            }
 
            // For every MeshPart specified in the list of slip boundaries, create a SlipFilterAssembler and
            // insert it into the std::map
            for(auto& it : slip_list)
            {
              // Create empty assembler
              auto new_asm = std::make_shared<Assembly::SlipFilterAssembler<TrafoType>>(this->_trafo);
 
              // Add the MeshPart to the assembler if it is there. There are legitimate reasons for it NOT to be
              // there, i.e. we are in parallel and our patch is not adjacent to that MeshPart
              auto* mpp = rmn_->find_mesh_part(it);
              if(mpp != nullptr)
              {
                new_asm->add_mesh_part(*mpp);
              }
 
              // Insert into the map
              String identifier(it);
              _slip_asm.emplace(identifier, new_asm);
            }
 
            for(Index i(0); i < _mu.size(); ++i)
            {
              _sum_mu += _mu(i);
            }
            sync_scalars.emplace("_sum_mu",&_sum_mu);
 
            if(mesh_conc_ != nullptr)
            {
              _mesh_conc = mesh_conc_->create_empty_clone();
              _mesh_conc->set_mesh_node(rmn_);
              _mesh_conc->add_sync_vecs(sync_vecs);
            }
 
            // Perform one time scal computation
            _compute_scales_once();
          }
 
        HyperelasticityFunctional() = delete;
        HyperelasticityFunctional(const HyperelasticityFunctional&) = delete;
 
        virtual ~HyperelasticityFunctional()
        {
        }
 
        bool empty() const
        {
          return (trafo_space.get_num_dofs() == Index(0));
        }
 
        VectorTypeL create_vector_l() const
        {
          return VectorTypeL(trafo_space.get_num_dofs());
        }
 
        VectorTypeR create_vector_r() const
        {
          return VectorTypeR(trafo_space.get_num_dofs());
        }
 
        Index columns() const
        {
          return _columns;
        }
 
        Index rows() const
        {
          return _rows;
        }
 
        static String name()
        {
          return "HyperelasticityFunctional<"+MeshType::name()+">";
        }
 
        virtual String info() const
        {
          const Index pad_width(30);
 
          String msg = name() + String("\nScale computation").pad_back(pad_width, '.')
            + String(": ") + stringify(_scale_computation) + "\n";
 
          msg += _cell_functional->info() + "\n";
          if(_mesh_conc != nullptr)
          {
            msg += _mesh_conc->info() + "\n";
          }
 
          return msg;
        }
 
        virtual void add_to_vtk_exporter(Geometry::ExportVTK<typename Trafo_::MeshType>& exporter) const override
        {
          exporter.add_cell_scalar("h", this->_h.elements());
          exporter.add_cell_scalar("lambda", this->_lambda.elements());
          exporter.add_cell_scalar("mu", this->_mu.elements());
 
          DataType* fval_norm(new DataType[this->get_mesh()->get_num_entities(MeshType::shape_dim)]);
          DataType* fval_det(new DataType[this->get_mesh()->get_num_entities(MeshType::shape_dim)]);
          DataType* fval_rec_det(new DataType[this->get_mesh()->get_num_entities(MeshType::shape_dim)]);
 
          DataType fval(0);
          eval_fval_cellwise(fval, fval_norm, fval_det, fval_rec_det);
          exporter.add_cell_vector("fval", fval_norm, fval_det, fval_rec_det);
 
          if(_mesh_conc != nullptr)
          {
            exporter.add_vertex_scalar("dist", _mesh_conc->get_dist().elements());
            exporter.add_vertex_vector("grad_dist", _mesh_conc->get_grad_dist());
 
            if(this->_penalty_param > DataType(0))
            {
              DataType* constraint_at_cell = new DataType[this->get_mesh()->get_num_entities(MeshType::shape_dim)];
 
              this->_mesh_conc->compute_constraint(constraint_at_cell);
              exporter.add_cell_scalar("alignment_constraint", constraint_at_cell);
 
              delete[] constraint_at_cell;
            }
          }
 
          delete [] fval_norm;
          delete [] fval_det;
          delete [] fval_rec_det;
 
        }
 
        DataType get_penalty_param() const
        {
          return _penalty_param;
        }
 
        void set_penalty_param(const DataType penalty_param_)
        {
          XASSERTM(penalty_param_ >= DataType(0),"Penalty parameter must be >= 0!");
          XASSERTM(penalty_param_ <= Math::pow(Math::huge<DataType>(), DataType(0.25)),
          "Excessively large penalty parameter.");
 
          _penalty_param = penalty_param_;
        }
 
        DataType get_constraint()
        {
          return _alignment_constraint;
        }
 
        //void set_mu(std::vector<CoordType>& cells, const CoordType& weight)
        //{
        //  for(Index i(0); i < cells.size(); ++i)
        //  {
        //    _mu(cells(i), weight);
        //  }
 
        //  _sum_mu = CoordType(0);
        //  for(Index i(0); i < _mu.size(); ++i)
        //  {
        //    _sum_mu += _mu(i);
        //  }
        //}
 
        virtual void init() override
        {
          init_pre_sync();
          init_post_sync();
        }
 
        virtual void init_pre_sync()
        {
          // Write any potential changes to the mesh
          this->buffer_to_mesh();
 
          // Compute distances if necessary
          if((this->_scale_computation == ScaleComputation::current_concentration) ||
              (this->_scale_computation == ScaleComputation::iter_concentration) ||
              (_penalty_param > DataType(0)) )
              {
                if(_mesh_conc == nullptr)
                  XABORTM("Scale computation set to "+stringify(_scale_computation)+"and alignment penalty factor to "+
                    stringify_fp_sci(_penalty_param)+", but no concentration function was given");
 
                _mesh_conc->compute_dist();
              }
 
          // Compute desired element size distribution
          this->_compute_scales_init();
        }
 
        virtual void init_post_sync()
        {
          // Scale mu so that || mu ||_1 = 1
          // For this, _sum_mu needs to have been summed up over all patches in the synchronization phase.
          if(_sum_mu != DataType(1))
          {
            _mu.scale(_mu, DataType(1)/_sum_mu);
            _sum_mu = DataType(1);
          }
 
          // Scale lambda so that || lambda ||_1 = 1
          // For this, _sum_lambda needs to have been summed up over all patches in the synchronization phase.
          if(_sum_lambda != DataType(1))
          {
            _lambda.scale(_lambda, DataType(1)/_sum_lambda);
            _sum_lambda = DataType(1);
          }
 
          // Compute the new optimal scales. For this, lambda needs to be scaled correctly and _sum_det needs to be
          // summed up over all patches in the synchronization phase.
          RefCellTrafo_::compute_h(_h, _lambda, _sum_det);
        }
 
        virtual void prepare(const VectorTypeR& vec_state, FilterType& filter)
        {
          prepare_pre_sync(vec_state, filter);
          prepare_post_sync(vec_state, filter);
        }
 
        virtual void prepare_pre_sync(const VectorTypeR& vec_state, FilterType& filter)
        {
          // Download state if necessary
          //CoordsBufferType vec_buf;
          //vec_buf.convert(vec_state);
 
          // Copy to buffer and mesh
          this->_coords_buffer.copy(vec_state);
          this->buffer_to_mesh();
 
          auto& dirichlet_filters = filter.template at<1>();
 
          for(auto& it : dirichlet_filters)
          {
            const auto& assembler = _dirichlet_asm.find(it.first);
 
            if(assembler == _dirichlet_asm.end())
              XABORTM("Could not find unit filter assembler for filter with key "+it.first);
 
            assembler->second->assemble(it.second, trafo_space, vec_state);
          }
 
          // The slip filter contains the outer unit normal, so reassemble it
          auto& slip_filters = filter.template at<0>();
          for(const auto& it:slip_filters)
            this->_mesh_node->adapt_by_name(it.first);
 
          // Copy back to the buffer as adaption might have changed the mesh
          this->mesh_to_buffer();
 
          for(auto& it : slip_filters)
          {
            const auto& assembler = _slip_asm.find(it.first);
 
            if(assembler == _slip_asm.end())
            {
              XABORTM("Could not find slip filter assembler for filter with key "+it.first);
            }
 
            assembler->second->assemble(it.second, trafo_space);
          }
 
          if( (this->_scale_computation == ScaleComputation::iter_concentration) ||
              (_penalty_param > DataType(0)))
              {
                if(_mesh_conc == nullptr)
                {
                  XABORTM("Scale computation set to "+stringify(_scale_computation)+"and alignment penalty factor to "+
                    stringify_fp_sci(_penalty_param)+", but no concentration function was given");
                }
 
                _mesh_conc->compute_dist();
              }
 
          _compute_scales_iter();
 
          if(_penalty_param > DataType(0))
          {
            _alignment_constraint = this->_mesh_conc->compute_constraint();
            this->sync_scalars.emplace("_alignment_constraint",&_alignment_constraint);
          }
 
        }
 
        virtual void prepare_post_sync(const VectorTypeR& DOXY(vec_state), FilterType& DOXY(filter))
        {
          if(_sum_lambda != DataType(1))
          {
            _lambda.scale(_lambda, DataType(1)/_sum_lambda);
            _sum_lambda = DataType(1);
          }
 
          RefCellTrafo_::compute_h(_h, _lambda, _sum_det);
 
          if(this->_scale_computation == ScaleComputation::iter_concentration)
          {
            _mesh_conc->compute_grad_h(this->_sum_det);
          }
        }
 
        virtual CoordType compute_cell_size_defect(CoordType& lambda_min, CoordType& lambda_max, CoordType& vol_min,
        CoordType& vol_max, CoordType& vol) const override
        {
 
          compute_cell_size_defect_pre_sync(vol_min, vol_max, vol);
          return compute_cell_size_defect_post_sync(lambda_min, lambda_max, vol_min, vol_max, vol);
 
        }
 
        void compute_cell_size_defect_pre_sync(CoordType& vol_min, CoordType& vol_max, CoordType& vol) const
        {
          vol = CoordType(0);
 
          vol_min = Math::huge<CoordType>();
          vol_max = CoordType(0);
 
          for(Index cell(0); cell < this->get_mesh()->get_num_entities(ShapeType::dimension); ++cell)
          {
            CoordType my_vol = this->_trafo.template compute_vol<ShapeType, CoordType>(cell);
            vol_min = Math::min(vol_min, my_vol);
            vol_max = Math::max(vol_min, my_vol);
            vol += my_vol;
          }
        }
 
        virtual CoordType compute_cell_size_defect_post_sync(CoordType& lambda_min, CoordType& lambda_max,
        CoordType& vol_min, CoordType& vol_max, const CoordType& vol) const
        {
          CoordType size_defect(0);
          lambda_min = Math::huge<CoordType>();
          lambda_max = CoordType(0);
 
          for(Index cell(0); cell < this->get_mesh()->get_num_entities(ShapeType::dimension); ++cell)
          {
            size_defect += Math::abs(this->_trafo.template
                compute_vol<ShapeType, CoordType>(cell)/vol - this->_lambda(cell));
 
            lambda_min = Math::min(lambda_min, this->_lambda(cell));
            lambda_max = Math::max(lambda_max, this->_lambda(cell));
          }
 
          vol_min /= vol;
          vol_max/= vol;
 
          return size_defect;
        }
 
        virtual void eval_fval_grad(CoordType& fval, VectorTypeL& grad, const bool& add_penalty_fval = true)
        {
          typedef typename TrafoType::template Evaluator<ShapeType, DataType>::Type TrafoEvaluator;
          typedef typename SpaceType::template Evaluator<TrafoEvaluator>::Type SpaceEvaluator;
 
          // Increase number of functional evaluations
          this->_num_func_evals++;
          this->_num_grad_evals++;
 
          // Total number of cells in the mesh
          Index ncells(this->get_mesh()->get_num_entities(ShapeType::dimension));
 
          // Index set for local/global numbering
          auto& idx = this->get_mesh()->template get_index_set<ShapeType::dimension,0>();
 
          // This will hold the coordinates for one element for passing to other routines
          FEAT::Tiny::Matrix <CoordType, Shape::FaceTraits<ShapeType,0>::count, MeshType::world_dim> x;
          // This will hold the local gradient for one element for passing to other routines
          FEAT::Tiny::Matrix<CoordType, Shape::FaceTraits<ShapeType,0>::count, MeshType::world_dim> grad_loc;
 
          // Clear gradient vector
          grad.format();
          fval = DataType(0);
 
          TrafoEvaluator trafo_eval(this->_trafo);
          SpaceEvaluator space_eval(this->trafo_space);
 
          // Compute the functional value for each cell
          for(Index cell(0); cell < ncells; ++cell)
          {
            DataType fval_loc(0);
            trafo_eval.prepare(cell);
            space_eval.prepare(trafo_eval);
 
            // Get local coordinates
            for(int j(0); j < Shape::FaceTraits<ShapeType,0>::count; ++j)
            {
              x[j] = this->_coords_buffer(idx(cell,j));
            }
 
            auto mat_tensor = RefCellTrafo_::compute_mat_tensor(x, this->_h(cell));
 
            this->_cell_functional->eval_fval_grad(
              fval_loc, grad_loc, mat_tensor, trafo_eval, space_eval, x, this->_h(cell));
 
            // Add the contribution from the dependence of h on the vertex coordinates
            if(this->_mesh_conc != nullptr && this->_mesh_conc->use_derivative())
            {
              const auto& grad_h = this->_mesh_conc->get_grad_h();
              this->_cell_functional->add_grad_h_part(
                grad_loc, mat_tensor, trafo_eval, space_eval, x, this->_h(cell), grad_h(cell));
            }
 
            fval += this->_mu(cell)*fval_loc;
            // Add local contributions to global gradient vector
            for(int j(0); j < Shape::FaceTraits<ShapeType,0>::count; ++j)
            {
              Index i(idx(cell,j));
              Tiny::Vector<CoordType, MeshType::world_dim, MeshType::world_dim> tmp(grad(i));
              tmp += this->_mu(cell)*grad_loc[j];
 
              grad(i,tmp);
            }
          }
 
          if(this->_penalty_param > DataType(0))
          {
            XASSERTM(this->_mesh_conc != nullptr,
            "You need a mesh concentration function for imposing a mesh alignment constraint!");
 
            if(add_penalty_fval)
            {
              // Add the quadratic penalty term to the functional value
              fval += this->_penalty_param*DataType(0.5)*Math::sqr(this->_alignment_constraint);
            }
 
            // Add the gradient of the penalty term
            this->_mesh_conc->add_constraint_grad(grad, this->_alignment_constraint, this->_penalty_param);
 
          }
        }
 
        virtual void eval_fval_cellwise(CoordType& fval, CoordType* fval_norm, CoordType* fval_cof,
        CoordType* fval_det) const
        {
          typedef typename TrafoType::template Evaluator<ShapeType, DataType>::Type TrafoEvaluator;
          typedef typename SpaceType::template Evaluator<TrafoEvaluator>::Type SpaceEvaluator;
 
          // Total number of cells in the mesh
          Index ncells(this->get_mesh()->get_num_entities(ShapeType::dimension));
 
          // Index set for local/global numbering
          auto& idx = this->get_mesh()->template get_index_set<ShapeType::dimension,0>();
 
          // This will hold the coordinates for one element for passing to other routines
          FEAT::Tiny::Matrix <CoordType, Shape::FaceTraits<ShapeType,0>::count, MeshType::world_dim> x;
          // Local cell dimensions for passing to other routines
          CoordType h(0);
          // This will hold the local gradient for one element for passing to other routines
          FEAT::Tiny::Matrix<CoordType, Shape::FaceTraits<ShapeType,0>::count, MeshType::world_dim> grad_loc;
 
          fval = DataType(0);
 
          TrafoEvaluator trafo_eval(this->_trafo);
          SpaceEvaluator space_eval(this->trafo_space);
 
          // Compute the functional value for each cell
          for(Index cell(0); cell < ncells; ++cell)
          {
            DataType fval_loc(0);
            trafo_eval.prepare(cell);
            space_eval.prepare(trafo_eval);
 
            h = this->_h(cell);
            // Get local coordinates
            for(int j(0); j < Shape::FaceTraits<ShapeType,0>::count; ++j)
            {
              x[j] = this->_coords_buffer(idx(cell,j));
            }
 
            auto mat_tensor = RefCellTrafo_::compute_mat_tensor(x, this->_h(cell));
 
            this->_cell_functional->eval_fval_cellwise(fval_loc, mat_tensor, trafo_eval, space_eval, x, h,
            fval_norm[cell], fval_cof[cell], fval_det[cell]);
 
            fval += this->_mu(cell)*fval_loc;
          }
 
        } // eval_fval_cellwise
 
      protected:
        virtual void _compute_lambda_cellsize()
        {
          Index ncells(this->get_mesh()->get_num_entities(ShapeType::dimension));
 
          _sum_lambda = DataType(0);
          for(Index cell(0); cell < ncells; ++cell)
          {
            _lambda(cell, this->_trafo.template compute_vol<ShapeType, CoordType>(cell));
            _sum_lambda += _lambda(cell);
          }
 
          this->sync_scalars.emplace("_sum_lambda",&_sum_lambda);
 
        }
 
        virtual void _compute_lambda_uniform()
        {
          _lambda.format(CoordType(1));
 
          _sum_lambda = CoordType(_lambda.size());
          this->sync_scalars.emplace("_sum_lambda",&_sum_lambda);
        }
 
        virtual void _compute_lambda_conc()
        {
          _mesh_conc->compute_conc();
          this->sync_scalars.emplace("_sum_conc",&(_mesh_conc->get_sum_conc()));
 
          _mesh_conc->compute_grad_conc();
 
          _mesh_conc->compute_grad_sum_det(this->_coords_buffer);
 
          _lambda.copy(_mesh_conc->get_conc());
 
          _sum_lambda = DataType(0);
          for(Index cell(0); cell < this->get_mesh()->get_num_entities(ShapeType::dimension); ++cell)
          {
            _sum_lambda += _lambda(cell);
          }
          this->sync_scalars.emplace("_sum_lambda",&_sum_lambda);
 
        }
 
        virtual void _compute_scales_init()
        {
 
          switch(_scale_computation)
          {
            case ScaleComputation::once_uniform:
            case ScaleComputation::once_cellsize:
            case ScaleComputation::once_concentration:
              return;
            case ScaleComputation::current_uniform:
              _compute_lambda_uniform();
              break;
            case ScaleComputation::current_cellsize:
              _compute_lambda_cellsize();
              break;
            case ScaleComputation::current_concentration:
            case ScaleComputation::iter_concentration:
              _compute_lambda_conc();
              break;
            default:
              return;
          }
 
          _sum_det = RefCellTrafo_::compute_sum_det(this->_coords_buffer, *(this->get_mesh()));
          this->sync_scalars.emplace("_sum_det",&_sum_det);
 
        }
 
        virtual void _compute_scales_iter()
        {
          switch(_scale_computation)
          {
            case ScaleComputation::iter_concentration:
              _compute_lambda_conc();
              break;
            default:
              return;
          }
 
          _sum_det = RefCellTrafo_::compute_sum_det(this->_coords_buffer, *(this->get_mesh()));
          this->sync_scalars.emplace("_sum_det",&_sum_det);
        }
 
        virtual void _compute_scales_once()
        {
          switch(_scale_computation)
          {
            case ScaleComputation::once_uniform:
              _compute_lambda_uniform();
              break;
            case ScaleComputation::once_cellsize:
              _compute_lambda_cellsize();
              break;
            case ScaleComputation::once_concentration:
              _compute_lambda_conc();
              break;
            default:
              return;
          }
 
          _sum_det = RefCellTrafo_::compute_sum_det(this->_coords_buffer, *(this->get_mesh()));
          this->sync_scalars.emplace("_sum_det",&_sum_det);
 
        }
 
    }; // class HyperelasticityFunctional
 
#ifdef FEAT_EICKT
    extern template class HyperelasticityFunctional
    <
      double, Index,
      Trafo::Standard::Mapping<Geometry::ConformalMesh< Shape::Simplex<2>, 2, double >>,
      Meshopt::RumpfFunctional
      <
        double,
        Trafo::Standard::Mapping<Geometry::ConformalMesh< Shape::Simplex<2>, 2, double >>
      >
    >;
 
    extern template class HyperelasticityFunctional
    <
      double, Index,
      Trafo::Standard::Mapping<Geometry::ConformalMesh< Shape::Simplex<3>, 3, double >>,
      Meshopt::RumpfFunctional
      <
        double,
        Trafo::Standard::Mapping<Geometry::ConformalMesh< Shape::Simplex<3>, 3, double >>
      >
    >;
 
    extern template class HyperelasticityFunctional
    <
      double, Index,
      Trafo::Standard::Mapping<Geometry::ConformalMesh< Shape::Hypercube<2>, 2, double >>,
      Meshopt::RumpfFunctional
      <
        double,
        Trafo::Standard::Mapping<Geometry::ConformalMesh< Shape::Hypercube<2>, 2, double >>
      >
    >;
 
    extern template class HyperelasticityFunctional
    <
      double, Index,
      Trafo::Standard::Mapping<Geometry::ConformalMesh< Shape::Hypercube<3>, 3, double >>,
      Meshopt::RumpfFunctional
      <
        double,
        Trafo::Standard::Mapping<Geometry::ConformalMesh< Shape::Hypercube<3>, 3, double >>
      >
    >;
#endif // FEAT_EICKT
 
  } // namespace Meshopt
} // namespace FEAT