Forked from
agnumpde / dune-tectonic
76 commits ahead of the upstream repository.
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
strikeslipfactory.cc 13.21 KiB
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/fufem/geometry/convexpolyhedron.hh>
#include <dune/contact/projections/normalprojection.hh>
#include "../problem-data/grid/cuboidgeometry.hh"
#include "../problem-data/bc.hh"
#include "../problem-data/myglobalfrictiondata.hh"
#include "../utils/diameter.hh"
#include "strikeslipfactory.hh"
template <class HostGridType, class VectorTEMPLATE>
void StrikeSlipFactory<HostGridType, VectorTEMPLATE>::setBodies() {
// set up cuboid geometries
std::array<double, 2> lengths = {this->parset_.template get<double>("body0.length"), this->parset_.template get<double>("body1.length")};
std::array<double, 2> heights = {this->parset_.template get<double>("body0.height"), this->parset_.template get<double>("body1.height")};
GlobalCoords origin(0);
const auto& frictionParset = this->parset_.sub("boundary.friction");
#if MY_DIM == 3 // TODO: not implemented
//std::array<double, 2> depths = {this->parset_.template get<double>("body0.depth"), this->parset_.template get<double>("body1.depth")};
//triangleGeometries_[0] = std::make_shared<TriangleGeometry>(origins[0], lengths[0], heights[0], depths[0]);
//triangleGeometries_[0]->addWeakeningPatch(frictionParset, {origins[0][0], origins[0][1]+ heights[0], 0}, {origins[0][0] + lengths[0], origins[0][1]+ heights[0], 0});
//triangleGeometries_[1] = std::make_shared<TriangleGeometry>(origins[1], lengths[1], heights[1], depths[1]);
//triangleGeometries_[1]->addWeakeningPatch(frictionParset, origins[1], {origins[1][0] + lengths[1], origins[1][1], 0});
#elif MY_DIM == 2
triangleGeometries_[0] = std::make_shared<TriangleGeometry>(origin, lengths[0], heights[0]);
triangleGeometries_[0]->addWeakeningPatch(frictionParset, triangleGeometries_[0]->A(), triangleGeometries_[0]->C());
triangleGeometries_[1] = std::make_shared<TriangleGeometry>(triangleGeometries_[0]->C(), -lengths[1], -heights[1]);
triangleGeometries_[1]->addWeakeningPatch(frictionParset, triangleGeometries_[1]->A(), triangleGeometries_[1]->C());
#else
#error CuboidGeometry only supports 2D and 3D!"
#endif
// set up reference grids
gridConstructor_ = std::make_unique<GridsConstructor<HostGridType>>(triangleGeometries_);
auto& grids = gridConstructor_->getGrids();
std::array<double, 2> smallestDiameter = {this->parset_.template get<double>("body0.smallestDiameter"), this->parset_.template get<double>("body1.smallestDiameter")};
for (size_t i=0; i<this->bodyCount_; i++) {
const auto& triangleGeometry = *triangleGeometries_[i];
// define weak patch and refine grid
const auto& weakeningRegions = triangleGeometry.weakeningRegions();
for (size_t j=0; j<weakeningRegions.size(); j++) {
refine(*grids[i], weakeningRegions[j], smallestDiameter[i], TriangleGeometry::lengthScale());
}
// determine minDiameter and maxDiameter
double minDiameter = std::numeric_limits<double>::infinity();
double maxDiameter = 0.0;
for (auto &&e : elements(grids[i]->leafGridView())) {
auto const geometry = e.geometry();
auto const diam = diameter(geometry);
minDiameter = std::min(minDiameter, diam);
maxDiameter = std::max(maxDiameter, diam);
}
std::cout << "Grid" << i << " min diameter: " << minDiameter << std::endl; std::cout << "Grid" << i << " max diameter: " << maxDiameter << std::endl;
}
#if MY_DIM == 2
bodyData_[0] = std::make_shared<MyBodyData<dim>>(this->parset_.sub("body0"), 0.0, zenith_());
this->bodies_[0] = std::make_shared<typename Base::LeafBody>(bodyData_[0], grids[0]);
bodyData_[1] = std::make_shared<MyBodyData<dim>>(this->parset_.sub("body1"), 0.0, zenith_());
this->bodies_[1] = std::make_shared<typename Base::LeafBody>(bodyData_[1], grids[1]);
#else
bodyData_[0] = std::make_shared<MyBodyData<dim>>(this->parset_.sub("body0"), this->parset_.template get<double>("gravity"), zenith_());
this->bodies_[0] = std::make_shared<typename Base::LeafBody>(bodyData_[0], grids[0]);
bodyData_[1] = std::make_shared<MyBodyData<dim>>(this->parset_.sub("body1"), this->parset_.template get<double>("gravity"), zenith_());
this->bodies_[1] = std::make_shared<typename Base::LeafBody>(bodyData_[1], grids[1]);
#endif
}
template <class HostGridType, class VectorTEMPLATE>
void StrikeSlipFactory<HostGridType, VectorTEMPLATE>::setLevelBodies() {
const size_t maxLevel = std::max({this->bodies_[0]->grid()->maxLevel(), this->bodies_[1]->grid()->maxLevel()});
for (size_t l=0; l<=maxLevel; l++) {
std::vector<size_t> bodyLevels(2, l);
this->contactNetwork_.addLevel(bodyLevels, l);
}
}
template <class HostGridType, class VectorTEMPLATE>
void StrikeSlipFactory<HostGridType, VectorTEMPLATE>::setCouplings() {
for (size_t i=0; i<this->bodyCount_; i++) {
const auto& triangleGeometry = *triangleGeometries_[i];
leafFaces_[i] = std::make_shared<LeafFaces>(this->bodies_[i]->gridView(), triangleGeometry);
levelFaces_[i] = std::make_shared<LevelFaces>(this->bodies_[i]->grid()->levelGridView(0), triangleGeometry);
}
auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<LeafBoundaryPatch>>();
nonmortarPatch_[0] = std::make_shared<LevelBoundaryPatch>(levelFaces_[1]->b);
mortarPatch_[0] = std::make_shared<LevelBoundaryPatch>(levelFaces_[0]->b);
weakPatches_[0] = std::make_shared<WeakeningRegion>(triangleGeometries_[1]->weakeningRegions()[0]);
auto& coupling = this->couplings_[0];
coupling = std::make_shared<typename Base::FrictionCouplingPair>();
coupling->set(1, 0, nonmortarPatch_[0], mortarPatch_[0], 0.1, Base::FrictionCouplingPair::CouplingType::STICK_SLIP, contactProjection, nullptr);
coupling->setWeakeningPatch(weakPatches_[0]);
coupling->setFrictionData(std::make_shared<MyGlobalFrictionData<GlobalCoords>>(this->parset_.sub("boundary.friction"), *weakPatches_[0], TriangleGeometry::lengthScale()));
}
template <class HostGridType, class VectorTEMPLATE>
void StrikeSlipFactory<HostGridType, VectorTEMPLATE>::setBoundaryConditions() {
using LeafBoundaryCondition = BoundaryCondition<LeafGridView, dim>;
using Function = Dune::VirtualFunction<double, double>;
std::shared_ptr<Function> neumannFunction = std::make_shared<NeumannCondition>();
std::shared_ptr<Function> velocityDirichletFunction = std::make_shared<VelocityDirichletCondition>(-1.0*this->parset_.template get<double>("boundary.dirichlet.finalVelocity"));
const double lengthScale = TriangleGeometry::lengthScale();
// body0
const auto& body0 = this->contactNetwork_.body(0);
const auto& leafVertexCount0 = body0->nVertices();
std::shared_ptr<Dune::BitSetVector<dim>> zeroDirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount0);
std::shared_ptr<Dune::BitSetVector<dim>> velocityDirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount0);
const auto& gridView0 = body0->gridView();
const auto& indexSet0 = gridView0.indexSet();
for (const auto& e : elements(gridView0)) { for (const auto& is : intersections(gridView0, e)) {
if (!is.boundary())
continue;
auto isCenter = is.geometry().center();
// Dirichlet boundary (c)
if (leafFaces_[0]->c.contains(is)) {
isCenter -= triangleGeometries_[0]->A();
if (isCenter.two_norm()/std::abs(triangleGeometries_[0]->length()) > 0.1) {
const auto inside = is.inside();
auto refElement = Dune::ReferenceElements<double, dim>::general(inside.type());
int n = refElement.size(is.indexInInside(),1,dim);
for (int i=0; i<n; i++) {
int faceIdxi = refElement.subEntity(is.indexInInside(), 1, i, dim);
int globalIdx = indexSet0.template subIndex(inside,faceIdxi,dim);
(*zeroDirichletNodes)[globalIdx][1] = true;
}
}
}
// velocity Dirichlet boundary (a)
if (leafFaces_[0]->a.contains(is)) {
isCenter -= triangleGeometries_[0]->C();
if (isCenter.two_norm()/std::abs(triangleGeometries_[0]->height()) > 0.1) {
const auto inside = is.inside();
auto refElement = Dune::ReferenceElements<double, dim>::general(inside.type());
int n = refElement.size(is.indexInInside(),1,dim);
for (int i=0; i<n; i++) {
int faceIdxi = refElement.subEntity(is.indexInInside(), 1, i, dim);
int globalIdx = indexSet0.template subIndex(inside,faceIdxi,dim);
(*velocityDirichletNodes)[globalIdx][0] = true;
}
}
}
}
}
std::shared_ptr<Function> zeroFunction = std::make_shared<NeumannCondition>();
std::shared_ptr<LeafBoundaryCondition> zeroDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
zeroDirichletBoundary->setBoundaryPatch(body0->gridView(), zeroDirichletNodes);
zeroDirichletBoundary->setBoundaryFunction(zeroFunction);
body0->addBoundaryCondition(zeroDirichletBoundary);
std::shared_ptr<LeafBoundaryCondition> velocityDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
velocityDirichletBoundary->setBoundaryPatch(body0->gridView(), velocityDirichletNodes);
velocityDirichletBoundary->setBoundaryFunction(velocityDirichletFunction);
body0->addBoundaryCondition(velocityDirichletBoundary);
// body1
const auto& body1 = this->contactNetwork_.body(1);
const auto& leafVertexCount1 = body1->nVertices();
std::shared_ptr<Dune::BitSetVector<dim>> zeroDirichletNodes1 = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount1);
std::shared_ptr<Dune::BitSetVector<dim>> loadNeumannNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount1);
const auto& gridView1 = body1->gridView();
const auto& indexSet1 = gridView1.indexSet();
for (const auto& e : elements(gridView1)) {
for (const auto& is : intersections(gridView1, e)) {
if (!is.boundary())
continue;
auto isCenter = is.geometry().center();
// Neumann boundary, normal load (c)
if (leafFaces_[1]->c.contains(is)) {
isCenter -= triangleGeometries_[1]->A();
if (isCenter.two_norm()/std::abs(triangleGeometries_[1]->length()) > 0.1) {
const auto inside = is.inside();
auto refElement = Dune::ReferenceElements<double, dim>::general(inside.type());
int n = refElement.size(is.indexInInside(),1,dim);
for (int i=0; i<n; i++) {
int faceIdxi = refElement.subEntity(is.indexInInside(), 1, i, dim);
int globalIdx = indexSet1.template subIndex(inside,faceIdxi,dim);
(*loadNeumannNodes)[globalIdx][1] = true;
}
}
}
// Dirichlet boundary (a)
if (leafFaces_[1]->a.contains(is)) {
isCenter -= triangleGeometries_[1]->C();
if (isCenter.two_norm()/std::abs(triangleGeometries_[1]->height()) > 0.1) {
const auto inside = is.inside();
auto refElement = Dune::ReferenceElements<double, dim>::general(inside.type());
int n = refElement.size(is.indexInInside(),1,dim);
for (int i=0; i<n; i++) {
int faceIdxi = refElement.subEntity(is.indexInInside(), 1, i, dim);
int globalIdx = indexSet1.template subIndex(inside,faceIdxi,dim);
(*zeroDirichletNodes1)[globalIdx][0] = true;
}
}
}
}
}
std::shared_ptr<Function> constantFunction = std::make_shared<NeumannCondition>(this->parset_.template get<double>("boundary.neumann.sigmaN"));
auto loadNeumannBoundary = std::make_shared<LeafBoundaryCondition>("neumann");
loadNeumannBoundary->setBoundaryPatch(body1->gridView(), loadNeumannNodes);
loadNeumannBoundary->setBoundaryFunction(constantFunction);
body1->addBoundaryCondition(loadNeumannBoundary);
std::shared_ptr<LeafBoundaryCondition> zeroDirichletBoundary1 = std::make_shared<LeafBoundaryCondition>("dirichlet");
zeroDirichletBoundary1->setBoundaryPatch(body1->gridView(), zeroDirichletNodes1);
zeroDirichletBoundary1->setBoundaryFunction(zeroFunction);
body1->addBoundaryCondition(zeroDirichletBoundary1);
// body0, body1: natural boundary conditions
for (size_t i=0; i<this->bodyCount_; i++) {
const auto& body = this->contactNetwork_.body(i);
// Neumann boundary
auto neumannBoundary = std::make_shared<LeafBoundaryCondition>("neumann");
auto neumannPatch = std::make_shared<LeafBoundaryPatch>(body->gridView(), true);
neumannBoundary->setBoundaryPatch(neumannPatch);
neumannBoundary->setBoundaryFunction(neumannFunction);
body->addBoundaryCondition(neumannBoundary);
}
}
#include "strikeslipfactory_tmpl.cc"