Commit e43ba4df authored by Sander, Oliver's avatar Sander, Oliver
Browse files

Use a composite basis for the brittle fracture simulation

Conceptually the brittle fracture simulation uses two spaces:
one for the deformation and one for the damage fields.
Previously, these were squeezed into a single power basis
with dim+1 coefficients.  This works, but it sometimes requires
manual index manipulations, which are not always obvious.
This patch introduces a proper composite basis that separates
the deformation from the damage.

The multi-indices for the global degrees of freedom remain unchanged!
parent e6b4dc05
Pipeline #4049 failed with stage
in 10 minutes and 32 seconds
......@@ -19,7 +19,7 @@ namespace Dune::FracturePhasefields
static constexpr std::size_t dim = Basis::GridView::dimension;
using Element = typename Basis::GridView::template Codim<0>::Entity;
using Tree = typename Basis::LocalView::Tree;
using FiniteElement = typename Tree::ChildType::FiniteElement;
using FiniteElement = typename std::tuple_element_t<0,typename Tree::ChildTypes>::ChildType::FiniteElement;
using DomainType = typename FiniteElement::Traits::LocalBasisType::Traits::DomainType;
using JacobianType = typename FiniteElement::Traits::LocalBasisType::Traits::JacobianType;
using JacobianFieldType = typename JacobianType::field_type;
......@@ -70,11 +70,10 @@ namespace Dune::FracturePhasefields
std::vector<FieldVector<double, dim*dim+1> >& jet,
const DomainType& position) const
{
const auto& fe = tree.child(0).finiteElement();
const auto& fe = tree.child(Indices::_0,0).finiteElement();
const auto& geometry = element.geometry();
const auto invJacobian = geometry.jacobianInverseTransposed(position);
using JacobianType = typename FiniteElement::Traits::LocalBasisType::Traits::JacobianType;
using RangeType = typename FiniteElement::Traits::LocalBasisType::Traits::RangeType;
auto values = std::vector<RangeType>(fe.localBasis().size());
......@@ -92,7 +91,7 @@ namespace Dune::FracturePhasefields
for(std::size_t i=0; i<dim; ++i)
{
auto deformationGradient = gradientTensor(gradient, i);
std::size_t index = tree.child(i).localIndex(j);
std::size_t index = tree.child(Indices::_0,i).localIndex(j);
for (std::size_t ii=0; ii<dim; ii++)
for (std::size_t jj=0; jj<dim; jj++)
jet[index][ii*dim+jj] = deformationGradient[ii][jj];
......@@ -100,7 +99,7 @@ namespace Dune::FracturePhasefields
}
// Evaluate shape function value
std::size_t index = tree.child(dim).localIndex(j);
std::size_t index = tree.child(Indices::_1,0).localIndex(j);
jet[index] = 0;
jet[index][dim*dim] = values[j];
}
......
......@@ -18,6 +18,7 @@
#include <dune/istl/bvector.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/functions/functionspacebases/compositebasis.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
......@@ -74,6 +75,7 @@
#include <dune/fracture-phasefields/damagedelasticenergydensity_nondecomposed.hh>
#include <dune/fracture-phasefields/damagedelasticenergydensity_largstrainstvenantkirchhoff.hh>
#include <dune/fracture-phasefields/damagedelasticenergydensity_wrapper.hh>
#include <dune/fracture-phasefields/transformedindexbasis.hh>
//#include <dune/matrix-vector/densematrixview.hh>
......@@ -353,28 +355,40 @@ int main (int argc, char *argv[]) try
using namespace Dune::Indices;
using namespace Dune::Functions::BasisBuilder;
// Create an index transformation that transforms the raw blocked indices
// to implement leaf blocking for composite nodes
auto transformation = Experimental::indexTransformation(
[](auto& multiIndex, const auto& basis) {
if (multiIndex[0] == 0)
multiIndex = {multiIndex[1], multiIndex[2]};
else if (multiIndex[0] == 1)
multiIndex = {multiIndex[1], dim};
},
[](const auto& prefix, const auto& basis) -> std::size_t {
if (prefix.size()==0)
return basis.size({0});
if (prefix.size()==1)
return dim+1;
return 0;
},
_1, _2); // Min/max index size
// A Dune::Functions basis spanning the whole ansatz space
auto basis = makeBasis(
grid->leafGridView(),
power<N>(
lagrange<1>(),
blockedInterleaved()
)
);
// Subspace bases for deformation and damage ansatz spaces
//
// Once dune-functions supports appropriate blocking,
// we can use
//
// Functions::subspaceBasis(basis, _0)
// Functions::subspaceBasis(basis, _0, i)
//
// to get bases of the full deformation space and the i-th
// component.
auto damageBasis = Functions::subspaceBasis(basis, Dune::index_constant<dim>());
// Vector marking all degerees of freedom used as Dirichlet nodes
Experimental::transformIndices(
composite(
power<dim>(lagrange<1>()),
lagrange<1>()
),
transformation));
using Basis = decltype(basis);
auto deformationBasis = Functions::subspaceBasis(basis, _0);
auto damageBasis = Functions::subspaceBasis(basis, _1);
// Vector marking all degrees of freedom used as Dirichlet nodes
Dune::BlockVector<Dune::FieldVector<bool,dim+1>> dirichletNodes(gridView.size(dim));
// Compute dirichlet nodes by evaluating a function given in python
......@@ -669,12 +683,7 @@ int main (int argc, char *argv[]) try
auto boundaryDeformationInput = Python::make_function<Coordinate>(pyModule.get("boundary_deformation"));
auto boundaryDeformationInputWrap = [&boundaryDeformationInput, &i](const FieldVector<double,dim>& x)
{
return boundaryDeformationInput(x, i);
};
Functions::interpolate(basis, Dune::Fufem::istlVectorBackend(x), boundaryDeformationInputWrap, Dune::Fufem::istlVectorBackend(dirichletNodes));
Functions::interpolate(deformationBasis, Fufem::istlVectorBackend(x), boundaryDeformationInput, Fufem::istlVectorBackend(dirichletNodes));
Vector lowerObstacle(x.size());
Vector upperObstacle(x.size());
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment