assemblers.cc

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <dune/istl/scaledidmatrix.hh>

#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/functions/constantfunction.hh>
#include <dune/fufem/assemblers/localassemblers/boundarymassassembler.hh>
#include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
#include <dune/fufem/assemblers/localassemblers/massassembler.hh>
#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
#include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
#include <dune/fufem/assemblers/localassemblers/viscosityassembler.hh>

#include <dune/tectonic/globalruinanonlinearity.hh>

#include "assemblers.hh"

template <class GridView, int dimension>
MyAssembler<GridView, dimension>::MyAssembler(GridView const &_gridView)
    : vertexBasis(_gridView),
      gridView(_gridView),
      vertexAssembler(vertexBasis, vertexBasis) {}

template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleFrictionalBoundaryMass(
    BoundaryPatch<GridView> const &frictionalBoundary,
    ScalarMatrix &frictionalBoundaryMass) {
  BoundaryMassAssembler<Grid, BoundaryPatch<GridView>, LocalVertexBasis,
                        LocalVertexBasis, Dune::FieldMatrix<double, 1, 1>> const
  frictionalBoundaryMassAssembler(frictionalBoundary);
  vertexAssembler.assembleOperator(frictionalBoundaryMassAssembler,
                                   frictionalBoundaryMass);
}

template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleMass(double density, Matrix &M) {
  MassAssembler<Grid, LocalVertexBasis, LocalVertexBasis,
                Dune::ScaledIdentityMatrix<double, dimension>> const localMass;
  vertexAssembler.assembleOperator(localMass, M);
  M *= density;
}

template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleElasticity(double E, double nu,
                                                          Matrix &A) {
  StVenantKirchhoffAssembler<Grid, LocalVertexBasis, LocalVertexBasis> const
  localStiffness(E, nu);
  vertexAssembler.assembleOperator(localStiffness, A);
}

template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleViscosity(double shearViscosity,
                                                         double bulkViscosity,
                                                         Matrix &C) {
  ViscosityAssembler<Grid, LocalVertexBasis, LocalVertexBasis> const
  localViscosity(shearViscosity, bulkViscosity);
  vertexAssembler.assembleOperator(localViscosity, C);
}

template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleBodyForce(double gravity,
                                                         double density,
                                                         LocalVector zenith,
                                                         Vector &f) {
  LocalVector weightedGravitationalDirection = zenith;
  weightedGravitationalDirection *= density * gravity;
  weightedGravitationalDirection *= -1;

  ConstantFunction<LocalVector, LocalVector> const gravityFunction(
      weightedGravitationalDirection);
  L2FunctionalAssembler<Grid, LocalVertexBasis, LocalVector>
  gravityFunctionalAssembler(gravityFunction);
  vertexAssembler.assembleFunctional(gravityFunctionalAssembler, f);
}

template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleNeumann(
    BoundaryPatch<GridView> const &neumannBoundary, Vector &f,
    Dune::VirtualFunction<double, double> const &neumann, double relativeTime) {
  LocalVector localNeumann(0);
  neumann.evaluate(relativeTime, localNeumann[0]);
  ConstantFunction<LocalVector, LocalVector> const fNeumann(localNeumann);
  NeumannBoundaryAssembler<Grid, LocalVector> neumannBoundaryAssembler(
      fNeumann);
  vertexAssembler.assembleBoundaryFunctional(neumannBoundaryAssembler, f,
                                             neumannBoundary);
}

template <class GridView, int dimension>
auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
    BoundaryPatch<GridView> const &frictionalBoundary, FrictionData const &fd)
    -> std::shared_ptr<GlobalNonlinearity<Matrix, Vector>> {
  // Lump negative normal stress (kludge)
  ScalarVector weights;
  {
    ConstantFunction<LocalVector, typename ScalarVector::block_type> const
    constantOneFunction(1);
    NeumannBoundaryAssembler<Grid, typename ScalarVector::block_type>
    frictionalBoundaryAssembler(constantOneFunction);
    vertexAssembler.assembleBoundaryFunctional(frictionalBoundaryAssembler,
                                               weights, frictionalBoundary);
  }
  for (size_t i = 0; i < weights.size(); ++i)
    assert(weights[i] >= 0.0);

  Dune::BitSetVector<1> frictionalNodes;
  frictionalBoundary.getVertices(frictionalNodes);
  return std::make_shared<GlobalRuinaNonlinearity<Matrix, Vector>>(
      frictionalNodes, weights, fd);
}

#include "assemblers_tmpl.cc"