Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • podlesny/dune-tectonic
  • agnumpde/dune-tectonic
2 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 1249 additions and 81 deletions
dune/tectonic/factories/contactnetworkfactory.hh 0 → 100644
View file @ 89161286
#ifndef SRC_CONTACTNETWORKFACTORY_HH
#define SRC_CONTACTNETWORKFACTORY_HH
#include <dune/common/parametertree.hh>
#include "../data-structures/network/contactnetwork.hh"
template <class HostGridType, class VectorType>
class ContactNetworkFactory {
public:
using ContactNetwork = ContactNetwork<HostGridType, VectorType>;
protected:
using LeafBody = typename ContactNetwork::LeafBody;
using FrictionCouplingPair = typename ContactNetwork::FrictionCouplingPair;
const Dune::ParameterTree& parset_;
const size_t bodyCount_;
const size_t couplingCount_;
std::vector<std::shared_ptr<LeafBody>> bodies_;
std::vector<std::shared_ptr<FrictionCouplingPair>> couplings_;
ContactNetwork contactNetwork_;
private:
virtual void setBodies() = 0;
virtual void setLevelBodies() = 0;
virtual void setCouplings() = 0;
virtual void setLevelCouplings() = 0;
virtual void setBoundaryConditions() = 0;
public:
ContactNetworkFactory(const Dune::ParameterTree& parset, size_t bodyCount, size_t couplingCount) :
parset_(parset),
bodyCount_(bodyCount),
couplingCount_(couplingCount),
bodies_(bodyCount_),
couplings_(couplingCount_),
contactNetwork_(bodyCount_, couplingCount_) {}
void build() {
setBodies();
contactNetwork_.setBodies(bodies_);
setCouplings();
contactNetwork_.setCouplings(couplings_);
setLevelBodies();
setLevelCouplings();
setBoundaryConditions();
}
ContactNetwork& contactNetwork() {
return contactNetwork_;
}
};
#endif
dune/tectonic/factories/levelcontactnetworkfactory.hh 0 → 100644
View file @ 89161286
#ifndef SRC_CONTACTNETWORKFACTORY_HH
#define SRC_CONTACTNETWORKFACTORY_HH
#include <dune/common/parametertree.hh>
#include "../data-structures/levelcontactnetwork.hh"
template <class HostGridType>
class ContactNetworkFactory {
public:
using LevelContactNetwork = LevelContactNetwork<GridType, dims>;
protected:
using Body = typename LevelContactNetwork::Body;
using FrictionCouplingPair = typename LevelContactNetwork::FrictionCouplingPair;
const Dune::ParameterTree& parset_;
const size_t bodyCount_;
const size_t couplingCount_;
std::vector<std::shared_ptr<Body>> bodies_;
std::vector<std::shared_ptr<FrictionCouplingPair>> couplings_;
LevelContactNetwork levelContactNetwork_;
private:
virtual void setBodies() = 0;
virtual void setCouplings() = 0;
virtual void setBoundaryConditions() = 0;
public:
LevelContactNetworkFactory(const Dune::ParameterTree& parset, size_t bodyCount, size_t couplingCount) :
parset_(parset),
bodyCount_(bodyCount),
couplingCount_(couplingCount),
bodies_(bodyCount_),
couplings_(couplingCount_),
levelContactNetwork_(bodyCount_, couplingCount_) {}
void build() {
setBodies();
levelContactNetwork_.setBodies(bodies_);
setCouplings();
levelContactNetwork_.setCouplings(couplings_);
setBoundaryConditions();
}
LevelContactNetwork& levelContactNetwork() {
return levelContactNetwork_;
}
};
#endif
dune/tectonic/factories/stackedblocksfactory.cc 0 → 100644
View file @ 89161286
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/fufem/geometry/convexpolyhedron.hh>
#include <dune/contact/projections/normalprojection.hh>
#include "../problem-data/bc.hh"
#include "../problem-data/grid/cuboidgeometry.hh"
#include "../problem-data/myglobalfrictiondata.hh"
#include "../utils/diameter.hh"
#include "stackedblocksfactory.hh"
template <class HostGridType, class VectorType>
void StackedBlocksFactory<HostGridType, VectorType>::setBodies() {
// set up cuboid geometries
double const length = 1.00;
double const height = 0.27;
double const weakLength = 0.60;
std::vector<GlobalCoords> origins(this->bodyCount_);
const auto& subParset = this->parset_.sub("boundary.friction.weakening");
#if MY_DIM == 3
double const depth = 0.60;
auto weakBound = [&] (const GlobalCoords& origin) {
GlobalCoords h = {origin[0] + weakLength, origin[1], origin[2]};
return h;
};
origins[0] = {0,0,0};
GlobalCoords lowerWeakOrigin = {0.2, 0, 0};
GlobalCoords upperWeakOrigin = {0.2, height, 0};
cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], length, height, depth);
cuboidGeometries_[0]->addWeakeningPatch(subParset, upperWeakOrigin, weakBound(upperWeakOrigin));
for (size_t i=1; i<this->bodyCount_-1; i++) {
origins[i] = cuboidGeometries_[i-1]->upperLeft();
GlobalCoords lowerWeakOrigin_i = lowerWeakOrigin + origins[i];
GlobalCoords upperWeakOrigin_i = upperWeakOrigin + origins[i];
cuboidGeometries_[i] = std::make_shared<CuboidGeometry>(origins[i], length, height, depth);
cuboidGeometries_[i]->addWeakeningPatch(subParset, upperWeakOrigin_i, weakBound(upperWeakOrigin_i));
cuboidGeometries_[i]->addWeakeningPatch(subParset, lowerWeakOrigin_i, weakBound(lowerWeakOrigin_i));
}
const size_t idx = this->bodyCount_-1;
origins[idx] = cuboidGeometries_[idx-1]->upperLeft();
lowerWeakOrigin += origins[idx];
cuboidGeometries_[idx] = std::make_shared<CuboidGeometry>(origins[idx], length, height, depth);
cuboidGeometries_[idx]->addWeakeningPatch(subParset, lowerWeakOrigin, weakBound(lowerWeakOrigin));
#elif MY_DIM == 2
auto weakBound = [&] (const GlobalCoords& origin) {
GlobalCoords h = {origin[0] + weakLength, origin[1]};
return h;
};
origins[0] = {0,0};
GlobalCoords lowerWeakOrigin = {0.2, 0};
GlobalCoords upperWeakOrigin = {0.2, height};
cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], length, height);
cuboidGeometries_[0]->addWeakeningPatch(subParset, upperWeakOrigin, weakBound(upperWeakOrigin));
for (size_t i=1; i<this->bodyCount_-1; i++) {
origins[i] = cuboidGeometries_[i-1]->upperLeft();
auto height_i = height/3.0;
GlobalCoords lowerWeakOrigin_i = lowerWeakOrigin + origins[i];
GlobalCoords upperWeakOrigin_i = {upperWeakOrigin[0], height_i};
upperWeakOrigin_i += origins[i];
cuboidGeometries_[i] = std::make_shared<CuboidGeometry>(origins[i], length, height_i);
cuboidGeometries_[i]->addWeakeningPatch(subParset, upperWeakOrigin_i, weakBound(upperWeakOrigin_i));
cuboidGeometries_[i]->addWeakeningPatch(subParset, lowerWeakOrigin_i, weakBound(lowerWeakOrigin_i));
}
const size_t idx = this->bodyCount_-1;
origins[idx] = cuboidGeometries_[idx-1]->upperLeft();
lowerWeakOrigin += origins[idx];
cuboidGeometries_[idx] = std::make_shared<CuboidGeometry>(origins[idx], length, height);
cuboidGeometries_[idx]->addWeakeningPatch(subParset, lowerWeakOrigin, weakBound(lowerWeakOrigin));
#else
#error CuboidGeometry only supports 2D and 3D!"
#endif
// set up reference grids
gridConstructor_ = std::make_unique<GridsConstructor<HostGridType>>(cuboidGeometries_);
auto& grids = gridConstructor_->getGrids();
for (size_t i=0; i<this->bodyCount_; i++) {
const auto& cuboidGeometry = *cuboidGeometries_[i];
// define weak patch and refine grid
const auto& weakeningRegions = cuboidGeometry.weakeningRegions();
for (size_t j=0; j<weakeningRegions.size(); j++) {
refine(*grids[i], weakeningRegions[j], this->parset_.template get<double>("boundary.friction.smallestDiameter"), CuboidGeometry::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;
}
for (size_t i=0; i<this->bodyCount_; i++) {
this->bodies_[i] = std::make_shared<typename Base::LeafBody>(bodyData_, grids[i]);
}
}
template <class HostGridType, class VectorType>
void StackedBlocksFactory<HostGridType, VectorType>::setCouplings() {
for (size_t i=0; i<this->bodyCount_; i++) {
const auto& cuboidGeometry = *cuboidGeometries_[i];
leafFaces_[i] = std::make_shared<LeafFaces>(this->bodies_[i]->gridView(), cuboidGeometry);
levelFaces_[i] = std::make_shared<LevelFaces>(this->bodies_[i]->grid()->levelGridView(0), cuboidGeometry);
}
auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<LeafBoundaryPatch>>();
std::shared_ptr<typename Base::FrictionCouplingPair::BackEndType> backend = nullptr;
for (size_t i=0; i<this->couplings_.size(); i++) {
auto& coupling = this->couplings_[i];
coupling = std::make_shared<typename Base::FrictionCouplingPair>();
nonmortarPatch_[i] = std::make_shared<LevelBoundaryPatch>(levelFaces_[i]->upper);
mortarPatch_[i] = std::make_shared<LevelBoundaryPatch>(levelFaces_[i+1]->lower);
weakPatches_[i] = std::make_shared<WeakeningRegion>(cuboidGeometries_[i]->weakeningRegions()[0]);
coupling->set(i, i+1, nonmortarPatch_[i], mortarPatch_[i], 0.1, Base::FrictionCouplingPair::CouplingType::STICK_SLIP, contactProjection, backend);
coupling->setWeakeningPatch(weakPatches_[i]);
coupling->setFrictionData(std::make_shared<MyGlobalFrictionData<GlobalCoords>>(this->parset_.sub("boundary.friction"), *weakPatches_[i], CuboidGeometry::lengthScale()));
}
}
template <class HostGridType, class VectorType>
void StackedBlocksFactory<HostGridType, VectorType>::setLevelBodies() {
const size_t nBodies = this->bodyCount_;
size_t maxLevel = 0;
std::vector<size_t> maxLevels(nBodies);
for (size_t i=0; i<nBodies; i++) {
const size_t bodyMaxLevel = this->bodies_[i]->grid()->maxLevel();
maxLevels[i] = bodyMaxLevel;
maxLevel = std::max(maxLevel, bodyMaxLevel);
}
for (int l=maxLevel; l>=0; l--) {
this->contactNetwork_.addLevel(maxLevels, l);
for (size_t i=0; i<nBodies; i++) {
maxLevels[i]--;
}
}
}
template <class HostGridType, class VectorType>
void StackedBlocksFactory<HostGridType, VectorType>::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>();
const double lengthScale = CuboidGeometry::lengthScale();
for (size_t i=0; i<this->bodyCount_; i++) {
const auto& body = this->contactNetwork_.body(i);
const auto& leafVertexCount = body->nVertices();
// Neumann boundary
std::shared_ptr<LeafBoundaryCondition> neumannBoundary = std::make_shared<LeafBoundaryCondition>(std::make_shared<LeafBoundaryPatch>(body->gridView()), neumannFunction, "neumann");
body->addBoundaryCondition(neumannBoundary);
// upper Dirichlet Boundary
// identify Dirichlet nodes on leaf level
if (i==this->bodyCount_-1) {
std::shared_ptr<Dune::BitSetVector<dim>> velocityDirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount);
for (int j=0; j<leafVertexCount; j++) {
if (leafFaces_[i]->upper.containsVertex(j))
(*velocityDirichletNodes)[j][0] = true;
}
std::shared_ptr<LeafBoundaryCondition> velocityDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
velocityDirichletBoundary->setBoundaryPatch(body->gridView(), velocityDirichletNodes);
velocityDirichletBoundary->setBoundaryFunction(velocityDirichletFunction);
body->addBoundaryCondition(velocityDirichletBoundary);
} else {
}
}
// lower Dirichlet Boundary
const auto& firstBody = this->contactNetwork_.body(0);
const auto& firstLeafVertexCount = firstBody->nVertices();
std::shared_ptr<Dune::BitSetVector<dim>> zeroDirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(firstLeafVertexCount);
for (int j=0; j<firstLeafVertexCount; j++) {
if (leafFaces_[0]->lower.containsVertex(j)) {
for (size_t d=0; d<dim; d++) {
(*zeroDirichletNodes)[j][d] = true;
}
}
#if MY_DIM == 3 //TODO: wrong, needs revision
if (leafFaces_[0]->front.containsVertex(j) || leafFaces_[0]->back.containsVertex(j))
(*zeroDirichletNodes)[j][2] = true;
#endif
}
std::shared_ptr<LeafBoundaryCondition> zeroDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
zeroDirichletBoundary->setBoundaryPatch(firstBody->gridView(), zeroDirichletNodes);
std::shared_ptr<Function> zeroFunction = std::make_shared<NeumannCondition>();
zeroDirichletBoundary->setBoundaryFunction(zeroFunction);
firstBody->addBoundaryCondition(zeroDirichletBoundary);
}
#include "stackedblocksfactory_tmpl.cc"
dune/tectonic/factories/stackedblocksfactory.hh 0 → 100644
View file @ 89161286
#ifndef SRC_STACKEDBLOCKSFACTORY_HH
#define SRC_STACKEDBLOCKSFACTORY_HH
#include <dune/common/bitsetvector.hh>
#include <dune/common/function.hh>
#include <dune/common/fvector.hh>
#include <dune/fufem/boundarypatch.hh>
#include "contactnetworkfactory.hh"
#include "../problem-data/mybody.hh"
#include "../problem-data/grid/mygrids.hh"
#include "../problem-data/grid/cuboidgeometry.hh"
template <class HostGridType, class VectorType> class StackedBlocksFactory : public ContactNetworkFactory<HostGridType, VectorType>{
private:
using Base = ContactNetworkFactory<HostGridType, VectorType>;
public:
using ContactNetwork = typename Base::ContactNetwork;
private:
using GlobalCoords = typename ContactNetwork::LocalVector;
using LeafGridView = typename ContactNetwork::GridView;
using LevelGridView = typename ContactNetwork::GridType::LevelGridView;
using LevelBoundaryPatch = BoundaryPatch<LevelGridView>;
using LeafBoundaryPatch = BoundaryPatch<LeafGridView>;
using LeafFaces = MyFaces<LeafGridView>;
using LevelFaces = MyFaces<LevelGridView>;
using CuboidGeometry= CuboidGeometry<typename GlobalCoords::field_type>;
using WeakeningRegion = typename CuboidGeometry::WeakeningRegion;
static const int dim = ContactNetwork::dim;
const std::shared_ptr<MyBodyData<dim>> bodyData_; // material properties of bodies
std::unique_ptr<GridsConstructor<HostGridType>> gridConstructor_;
std::vector<std::shared_ptr<CuboidGeometry>> cuboidGeometries_;
std::vector<std::shared_ptr<LeafFaces>> leafFaces_;
std::vector<std::shared_ptr<LevelFaces>> levelFaces_;
std::vector<std::shared_ptr<WeakeningRegion>> weakPatches_;
std::vector<std::shared_ptr<LevelBoundaryPatch>> nonmortarPatch_;
std::vector<std::shared_ptr<LevelBoundaryPatch>> mortarPatch_;
public:
StackedBlocksFactory(const Dune::ParameterTree& parset) :
Base(parset, parset.get<size_t>("problem.bodyCount"), parset.get<size_t>("problem.bodyCount")-1),
bodyData_(std::make_shared<MyBodyData<dim>>(this->parset_.sub("body"), this->parset_.template get<double>("gravity"), zenith_())),
cuboidGeometries_(this->bodyCount_),
leafFaces_(this->bodyCount_),
levelFaces_(this->bodyCount_),
weakPatches_(this->bodyCount_),
nonmortarPatch_(this->couplingCount_),
mortarPatch_(this->couplingCount_)
{}
void setBodies();
void setLevelBodies();
void setCouplings();
void setLevelCouplings() {}
void setBoundaryConditions();
auto& weakPatches() {
return weakPatches_;
}
private:
static constexpr Dune::FieldVector<double, MY_DIM> zenith_() {
#if MY_DIM == 2
return {0, 1};
#elif MY_DIM == 3
return {0, 1, 0};
#endif
}
};
#endif
src/hdf5/surface-writer_tmpl.cc → dune/tectonic/factories/stackedblocksfactory_tmpl.cc
View file @ 89161286
#include "../explicitvectors.hh" #ifndef MY_DIM
#include "../explicitgrid.hh" #error MY_DIM unset
#endif
#include "../program_state.hh"
using MyProgramState = ProgramState<Vector, ScalarVector>; #include "../explicitgrid.hh"
#include "../explicitvectors.hh"
template class SurfaceWriter<MyProgramState, GridView>; template class StackedBlocksFactory<Grid, Vector>;
dune/tectonic/factories/threeblocksfactory.cc 0 → 100644
View file @ 89161286
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/fufem/geometry/convexpolyhedron.hh>
#include <dune/contact/projections/normalprojection.hh>
#include "../problem-data/bc.hh"
#include "../problem-data/myglobalfrictiondata.hh"
#include "../utils/diameter.hh"
#include "threeblocksfactory.hh"
template <class HostGridType, class VectorTEMPLATE>
void ThreeBlocksFactory<HostGridType, VectorTEMPLATE>::setBodies() {
// set up cuboid geometries
std::array<double, 3> lengths = {3.0, 1.0, 1.0};
std::array<double, 3> heights = {1.0, 1.0, 1.0};
std::array<GlobalCoords, 3> origins;
const auto& subParset = this->parset_.sub("boundary.friction.weakening");
#if MY_DIM == 3 // TODO: not implemented
double const depth = 0.60;
origins[0] = {0, 0, 0};
origins[1] = {0, origins[0][1] + heights[0], 0};
origins[2] = {origins[1][0] + lengths[1], origins[0][1] + heights[0], 0};
cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], lengths[0], heights[0], depth);
cuboidGeometries_[0]->addWeakeningPatch(subParset, {origins[0][0], origins[0][1]+ heights[0], 0}, {origins[0][0] + lengths[0], origins[0][1]+ heights[0], 0});
cuboidGeometries_[1] = std::make_shared<CuboidGeometry>(origins[1], lengths[1], heights[1], depth);
cuboidGeometries_[1]->addWeakeningPatch(subParset, origins[1], {origins[1][0] + lengths[1], origins[1][1], 0});
cuboidGeometries_[1]->addWeakeningPatch(subParset, origins[1], {origins[1][0] + lengths[1], origins[1][1] + heights[1], 0});
cuboidGeometries_[2] = std::make_shared<CuboidGeometry>(origins[2], lengths[2], heights[2], depth);
cuboidGeometries_[2]->addWeakeningPatch(subParset, origins[2], {origins[2][0] + lengths[2], origins[2][1], 0});
cuboidGeometries_[2]->addWeakeningPatch(subParset, origins[2], {origins[2][0], origins[2][1] + widths[2], 0});
#elif MY_DIM == 2
origins[0] = {0, 0};
origins[1] = {0, origins[0][1] + heights[0]};
origins[2] = {origins[1][0] + lengths[1], origins[0][1] + heights[0]};
cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], lengths[0], heights[0]);
cuboidGeometries_[0]->addWeakeningPatch(subParset, {origins[0][0], origins[0][1]+ heights[0]}, {origins[0][0] + lengths[0], origins[0][1]+ heights[0]});
cuboidGeometries_[1] = std::make_shared<CuboidGeometry>(origins[1], lengths[1], heights[1]);
cuboidGeometries_[1]->addWeakeningPatch(subParset, origins[1], {origins[1][0] + lengths[1], origins[1][1]});
cuboidGeometries_[1]->addWeakeningPatch(subParset, {origins[1][0] + lengths[1], origins[1][1]}, {origins[1][0] + lengths[1], origins[1][1] + heights[1]});
cuboidGeometries_[2] = std::make_shared<CuboidGeometry>(origins[2], lengths[2], heights[2]);
cuboidGeometries_[2]->addWeakeningPatch(subParset, origins[2], {origins[2][0] + lengths[2], origins[2][1]});
cuboidGeometries_[2]->addWeakeningPatch(subParset, origins[2], {origins[2][0], origins[2][1] + heights[2]});
#else
#error CuboidGeometry only supports 2D and 3D!"
#endif
// set up reference grids
gridConstructor_ = std::make_unique<GridsConstructor<HostGridType>>(cuboidGeometries_);
auto& grids = gridConstructor_->getGrids();
for (size_t i=0; i<this->bodyCount_; i++) {
const auto& cuboidGeometry = *cuboidGeometries_[i];
// define weak patch and refine grid
const auto& weakeningRegions = cuboidGeometry.weakeningRegions();
for (size_t j=0; j<weakeningRegions.size(); j++) {
refine(*grids[i], weakeningRegions[j], this->parset_.template get<double>("boundary.friction.smallestDiameter"), CuboidGeometry::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;
}
for (size_t i=0; i<this->bodyCount_; i++) {
this->bodies_[i] = std::make_shared<typename Base::LeafBody>(bodyData_, grids[i]);
}
}
template <class HostGridType, class VectorTEMPLATE>
void ThreeBlocksFactory<HostGridType, VectorTEMPLATE>::setLevelBodies() {
const size_t maxLevel = std::max({this->bodies_[0]->grid()->maxLevel(), this->bodies_[1]->grid()->maxLevel(), this->bodies_[2]->grid()->maxLevel()});
for (size_t l=0; l<=maxLevel; l++) {
std::vector<size_t> bodyLevels(3, l);
this->contactNetwork_.addLevel(bodyLevels, l);
}
}
template <class HostGridType, class VectorTEMPLATE>
void ThreeBlocksFactory<HostGridType, VectorTEMPLATE>::setCouplings() {
for (size_t i=0; i<this->bodyCount_; i++) {
const auto& cuboidGeometry = *cuboidGeometries_[i];
std::cout << this->bodies_[i]->nVertices() << std::endl;
std::cout << (leafFaces_[i] == nullptr) << std::endl;
leafFaces_[i] = std::make_shared<LeafFaces>(this->bodies_[i]->gridView(), cuboidGeometry);
levelFaces_[i] = std::make_shared<LevelFaces>(this->bodies_[i]->grid()->levelGridView(0), cuboidGeometry);
}
auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<LeafBoundaryPatch>>();
std::shared_ptr<typename Base::FrictionCouplingPair::BackEndType> backend = nullptr;
std::array<std::array<size_t, 2>, 3> couplingIndices;
nonmortarPatch_[0] = std::make_shared<LevelBoundaryPatch>(levelFaces_[1]->lower);
mortarPatch_[0] = std::make_shared<LevelBoundaryPatch>(levelFaces_[0]->upper);
weakPatches_[0] = std::make_shared<WeakeningRegion>(cuboidGeometries_[1]->weakeningRegions()[0]);
couplingIndices[0] = {1, 0};
nonmortarPatch_[1] = std::make_shared<LevelBoundaryPatch>(levelFaces_[2]->lower);
mortarPatch_[1] = std::make_shared<LevelBoundaryPatch>(levelFaces_[0]->upper);
weakPatches_[1] = std::make_shared<WeakeningRegion>(cuboidGeometries_[2]->weakeningRegions()[0]);
couplingIndices[1] = {2, 0};
nonmortarPatch_[2] = std::make_shared<LevelBoundaryPatch>(levelFaces_[1]->right);
mortarPatch_[2] = std::make_shared<LevelBoundaryPatch>(levelFaces_[2]->left);
weakPatches_[2] = std::make_shared<WeakeningRegion>(cuboidGeometries_[1]->weakeningRegions()[1]);
couplingIndices[2] = {1, 2};
for (size_t i=0; i<this->couplings_.size(); i++) {
auto& coupling = this->couplings_[i];
coupling = std::make_shared<typename Base::FrictionCouplingPair>();
coupling->set(couplingIndices[i][0], couplingIndices[i][1], nonmortarPatch_[i], mortarPatch_[i], 0.1, Base::FrictionCouplingPair::CouplingType::STICK_SLIP, contactProjection, backend);
coupling->setWeakeningPatch(weakPatches_[i]);
coupling->setFrictionData(std::make_shared<MyGlobalFrictionData<GlobalCoords>>(this->parset_.sub("boundary.friction"), *weakPatches_[i], CuboidGeometry::lengthScale()));
}
}
template <class HostGridType, class VectorTEMPLATE>
void ThreeBlocksFactory<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>();
const double lengthScale = CuboidGeometry::lengthScale();
for (size_t i=0; i<this->bodyCount_; i++) {
const auto& body = this->contactNetwork_.body(i);
const auto& leafVertexCount = body->nVertices();
// Neumann boundary
std::shared_ptr<LeafBoundaryCondition> neumannBoundary = std::make_shared<LeafBoundaryCondition>(std::make_shared<LeafBoundaryPatch>(body->gridView()), neumannFunction, "neumann");
body->addBoundaryCondition(neumannBoundary);
// upper Dirichlet Boundary
// identify Dirichlet nodes on leaf level
if (i>0) {
std::shared_ptr<Dune::BitSetVector<dim>> velocityDirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount);
for (int j=0; j<leafVertexCount; j++) {
if (leafFaces_[i]->upper.containsVertex(j))
(*velocityDirichletNodes)[j][0] = true;
#if MY_DIM == 3 //TODO: wrong, needs revision
if (leafFaces_[i]->front.containsVertex(j) || leafFaces_[i]->back.containsVertex(j))
zeroDirichletNodes->at(j)[2] = true;
#endif
}
std::shared_ptr<LeafBoundaryCondition> velocityDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
velocityDirichletBoundary->setBoundaryPatch(body->gridView(), velocityDirichletNodes);
velocityDirichletBoundary->setBoundaryFunction(velocityDirichletFunction);
body->addBoundaryCondition(velocityDirichletBoundary);
}
}
// lower Dirichlet Boundary
const auto& firstBody = this->contactNetwork_.body(0);
const auto& firstLeafVertexCount = firstBody->nVertices();
std::shared_ptr<Dune::BitSetVector<dim>> zeroDirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(firstLeafVertexCount);
for (int j=0; j<firstLeafVertexCount; j++) {
if (leafFaces_[0]->lower.containsVertex(j)) {
for (size_t d=0; d<dim; d++) {
(*zeroDirichletNodes)[j][d] = true;
}
}
}
std::shared_ptr<LeafBoundaryCondition> zeroDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
zeroDirichletBoundary->setBoundaryPatch(firstBody->gridView(), zeroDirichletNodes);
std::shared_ptr<Function> zeroFunction = std::make_shared<NeumannCondition>();
zeroDirichletBoundary->setBoundaryFunction(zeroFunction);
firstBody->addBoundaryCondition(zeroDirichletBoundary);
}
#include "threeblocksfactory_tmpl.cc"
dune/tectonic/factories/threeblocksfactory.hh 0 → 100644
View file @ 89161286
#ifndef SRC_THREEBLOCKSFACTORY_HH
#define SRC_THREEBLOCKSFACTORY_HH
#include <dune/common/bitsetvector.hh>
#include <dune/common/function.hh>
#include <dune/common/fvector.hh>
#include <dune/fufem/boundarypatch.hh>
#include "contactnetworkfactory.hh"
#include "../problem-data/mybody.hh"
#include "../problem-data/grid/mygrids.hh"
#include "../problem-data/grid/cuboidgeometry.hh"
template <class HostGridType, class VectorType>
class ThreeBlocksFactory : public ContactNetworkFactory<HostGridType, VectorType>{
private:
using Base = ContactNetworkFactory<HostGridType, VectorType>;
public:
using ContactNetwork = typename Base::ContactNetwork;
private:
using GlobalCoords = typename ContactNetwork::LocalVector;
using LeafGridView = typename ContactNetwork::GridView;
using LevelGridView = typename ContactNetwork::GridType::LevelGridView;
using LevelBoundaryPatch = BoundaryPatch<LevelGridView>;
using LeafBoundaryPatch = BoundaryPatch<LeafGridView>;
using LeafFaces = MyFaces<LeafGridView>;
using LevelFaces = MyFaces<LevelGridView>;
using CuboidGeometry= CuboidGeometry<typename GlobalCoords::field_type>;
using WeakeningRegion = typename CuboidGeometry::WeakeningRegion;
static const int dim = ContactNetwork::dim;
const std::shared_ptr<MyBodyData<dim>> bodyData_; // material properties of bodies
std::unique_ptr<GridsConstructor<HostGridType>> gridConstructor_;
std::vector<std::shared_ptr<CuboidGeometry>> cuboidGeometries_;
std::vector<std::shared_ptr<LeafFaces>> leafFaces_;
std::vector<std::shared_ptr<LevelFaces>> levelFaces_;
std::vector<std::shared_ptr<WeakeningRegion>> weakPatches_;
std::vector<std::shared_ptr<LevelBoundaryPatch>> nonmortarPatch_;
std::vector<std::shared_ptr<LevelBoundaryPatch>> mortarPatch_;
public:
ThreeBlocksFactory(const Dune::ParameterTree& parset) :
Base(parset, 3, 3),
bodyData_(std::make_shared<MyBodyData<dim>>(this->parset_.sub("body"), this->parset_.template get<double>("gravity"), zenith_())),
cuboidGeometries_(3),
leafFaces_(3),
levelFaces_(3),
weakPatches_(3),
nonmortarPatch_(3),
mortarPatch_(3)
{}
void setBodies();
void setLevelBodies();
void setCouplings();
void setLevelCouplings() {}
void setBoundaryConditions();
auto& weakPatches() {
return weakPatches_;
}
private:
static constexpr Dune::FieldVector<double, MY_DIM> zenith_() {
#if MY_DIM == 2
return {0, 1};
#elif MY_DIM == 3
return {0, 1, 0};
#endif
}
};
#endif
dune/tectonic/factories/threeblocksfactory_tmpl.cc 0 → 100644
View file @ 89161286
#ifndef MY_DIM
#error MY_DIM unset
#endif
#include "../explicitgrid.hh"
#include "../explicitvectors.hh"
template class ThreeBlocksFactory<Grid, Vector>;
dune/tectonic/factories/twoblocksfactory.cc 0 → 100644
View file @ 89161286
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/fufem/geometry/convexpolyhedron.hh>
#include <dune/contact/projections/normalprojection.hh>
#include "../problem-data/bc.hh"
#include "../problem-data/myglobalfrictiondata.hh"
#include "../utils/diameter.hh"
#include "twoblocksfactory.hh"
template <class HostGridType, class VectorTEMPLATE>
void TwoBlocksFactory<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")};
std::array<GlobalCoords, 2> origins;
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")};
origins[0] = {0, 0, 0};
origins[1] = {lengths[0]/2.0, origins[0][1] + heights[0], 0};
cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], lengths[0], heights[0], depths[0]);
cuboidGeometries_[0]->addWeakeningPatch(frictionParset, {origins[0][0], origins[0][1]+ heights[0], 0}, {origins[0][0] + lengths[0], origins[0][1]+ heights[0], 0});
cuboidGeometries_[1] = std::make_shared<CuboidGeometry>(origins[1], lengths[1], heights[1], depths[1]);
cuboidGeometries_[1]->addWeakeningPatch(frictionParset, origins[1], {origins[1][0] + lengths[1], origins[1][1], 0});
#elif MY_DIM == 2
origins[0] = {0, 0};
origins[1] = {lengths[0]/2.0, origins[0][1] + heights[0]};
cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], lengths[0], heights[0]);
cuboidGeometries_[0]->addWeakeningPatch(frictionParset, {origins[0][0], origins[0][1]+ heights[0]}, {origins[0][0] + lengths[0], origins[0][1]+ heights[0]});
cuboidGeometries_[1] = std::make_shared<CuboidGeometry>(origins[1], lengths[1], heights[1]);
cuboidGeometries_[1]->addWeakeningPatch(frictionParset, {origins[1][0] + 0.2*lengths[1], origins[1][1]}, {origins[1][0] + 0.8*lengths[1], origins[1][1]});
#else
#error CuboidGeometry only supports 2D and 3D!"
#endif
// set up reference grids
gridConstructor_ = std::make_unique<GridsConstructor<HostGridType>>(cuboidGeometries_);
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& cuboidGeometry = *cuboidGeometries_[i];
// define weak patch and refine grid
const auto& weakeningRegions = cuboidGeometry.weakeningRegions();
for (size_t j=0; j<weakeningRegions.size(); j++) {
refine(*grids[i], weakeningRegions[j], smallestDiameter[i], CuboidGeometry::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;
}
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]);
}
template <class HostGridType, class VectorTEMPLATE>
void TwoBlocksFactory<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 TwoBlocksFactory<HostGridType, VectorTEMPLATE>::setCouplings() {
for (size_t i=0; i<this->bodyCount_; i++) {
const auto& cuboidGeometry = *cuboidGeometries_[i];
leafFaces_[i] = std::make_shared<LeafFaces>(this->bodies_[i]->gridView(), cuboidGeometry);
levelFaces_[i] = std::make_shared<LevelFaces>(this->bodies_[i]->grid()->levelGridView(0), cuboidGeometry);
}
auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<LeafBoundaryPatch>>();
nonmortarPatch_[0] = std::make_shared<LevelBoundaryPatch>(levelFaces_[1]->lower);
mortarPatch_[0] = std::make_shared<LevelBoundaryPatch>(levelFaces_[0]->upper);
weakPatches_[0] = std::make_shared<WeakeningRegion>(cuboidGeometries_[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], CuboidGeometry::lengthScale()));
}
template <class HostGridType, class VectorTEMPLATE>
void TwoBlocksFactory<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>();
const double lengthScale = CuboidGeometry::lengthScale();
// body0: Dirichlet boundary (lower)
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);
for (int j=0; j<leafVertexCount0; j++) {
if (leafFaces_[0]->lower.containsVertex(j)) {
for (size_t d=0; d<dim; d++) {
(*zeroDirichletNodes)[j][d] = true;
}
}
}
std::shared_ptr<LeafBoundaryCondition> zeroDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
zeroDirichletBoundary->setBoundaryPatch(body0->gridView(), zeroDirichletNodes);
std::shared_ptr<Function> zeroFunction = std::make_shared<NeumannCondition>();
zeroDirichletBoundary->setBoundaryFunction(zeroFunction);
body0->addBoundaryCondition(zeroDirichletBoundary);
// body1: Dirichlet boundary (upper)
const auto& body1 = this->contactNetwork_.body(1);
const auto& leafVertexCount1 = body1->nVertices();
std::shared_ptr<Dune::BitSetVector<dim>> velocityDirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount1);
for (int j=0; j<leafVertexCount1; j++) {
if (leafFaces_[1]->upper.containsVertex(j))
(*velocityDirichletNodes)[j][0] = true;
#if MY_DIM == 3 //TODO: wrong, needs revision
if (leafFaces_[1]->front.containsVertex(j) || leafFaces_[1]->back.containsVertex(j))
zeroDirichletNodes->at(j)[2] = true;
#endif
}
std::shared_ptr<LeafBoundaryCondition> velocityDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
velocityDirichletBoundary->setBoundaryPatch(body1->gridView(), velocityDirichletNodes);
velocityDirichletBoundary->setBoundaryFunction(velocityDirichletFunction);
body1->addBoundaryCondition(velocityDirichletBoundary);
// body0, body1: natural boundary conditions
for (size_t i=0; i<this->bodyCount_; i++) {
const auto& body = this->contactNetwork_.body(i);
const auto& leafVertexCount = body->nVertices();
// Neumann boundary
std::shared_ptr<LeafBoundaryCondition> neumannBoundary = std::make_shared<LeafBoundaryCondition>(std::make_shared<LeafBoundaryPatch>(body->gridView()), neumannFunction, "neumann");
body->addBoundaryCondition(neumannBoundary);
}
// body1: Neumann boundary (upper)
// normal load
/*std::shared_ptr<Dune::BitSetVector<dim>> loadNeumannNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount1);
for (int j=0; j<leafVertexCount1; j++) {
if (leafFaces_[1]->upper.containsVertex(j))
(*loadNeumannNodes)[j][1] = true;
#if MY_DIM == 3 //TODO: wrong, needs revision
if (leafFaces_[1]->front.containsVertex(j) || leafFaces_[1]->back.containsVertex(j))
zeroDirichletNodes->at(j)[2] = true;
#endif
}
std::shared_ptr<LeafBoundaryCondition> loadNeumannBoundary = std::make_shared<LeafBoundaryCondition>("neumann");
std::shared_ptr<Function> constantFunction = std::make_shared<NeumannCondition>(this->parset_.template get<double>("boundary.neumann.sigmaN"));
loadNeumannBoundary->setBoundaryPatch(body1->gridView(), loadNeumannNodes);
loadNeumannBoundary->setBoundaryFunction(constantFunction);
body1->addBoundaryCondition(loadNeumannBoundary);*/
}
#include "twoblocksfactory_tmpl.cc"
dune/tectonic/factories/twoblocksfactory.hh 0 → 100644
View file @ 89161286
#ifndef DUNE_TECTONIC_FACTORIES_TWOBLOCKSFACTORY_HH
#define DUNE_TECTONIC_FACTORIES_TWOBLOCKSFACTORY_HH
#include <dune/common/bitsetvector.hh>
#include <dune/common/function.hh>
#include <dune/common/fvector.hh>
#include <dune/fufem/boundarypatch.hh>
#include "contactnetworkfactory.hh"
#include "../problem-data/mybody.hh"
#include "../problem-data/grid/mygrids.hh"
#include "../problem-data/grid/cuboidgeometry.hh"
template <class HostGridType, class VectorType> class TwoBlocksFactory : public ContactNetworkFactory<HostGridType, VectorType>{
private:
using Base = ContactNetworkFactory<HostGridType, VectorType>;
public:
using ContactNetwork = typename Base::ContactNetwork;
private:
using GlobalCoords = typename ContactNetwork::LocalVector;
using LeafGridView = typename ContactNetwork::GridView;
using LevelGridView = typename ContactNetwork::GridType::LevelGridView;
using LevelBoundaryPatch = BoundaryPatch<LevelGridView>;
using LeafBoundaryPatch = BoundaryPatch<LeafGridView>;
using LeafFaces = MyFaces<LeafGridView>;
using LevelFaces = MyFaces<LevelGridView>;
using CuboidGeometry= CuboidGeometry<typename GlobalCoords::field_type>;
using WeakeningRegion = typename CuboidGeometry::WeakeningRegion;
static const int dim = ContactNetwork::dim;
std::vector<std::shared_ptr<MyBodyData<dim>>> bodyData_; // material properties of bodies
std::unique_ptr<GridsConstructor<HostGridType>> gridConstructor_;
std::vector<std::shared_ptr<CuboidGeometry>> cuboidGeometries_;
std::vector<std::shared_ptr<LeafFaces>> leafFaces_;
std::vector<std::shared_ptr<LevelFaces>> levelFaces_;
std::vector<std::shared_ptr<WeakeningRegion>> weakPatches_;
std::vector<std::shared_ptr<LevelBoundaryPatch>> nonmortarPatch_;
std::vector<std::shared_ptr<LevelBoundaryPatch>> mortarPatch_;
public:
TwoBlocksFactory(const Dune::ParameterTree& parset) :
Base(parset, 2, 1),
bodyData_(2),
cuboidGeometries_(2),
leafFaces_(2),
levelFaces_(2),
weakPatches_(2),
nonmortarPatch_(1),
mortarPatch_(1)
{}
void setBodies();
void setLevelBodies();
void setCouplings();
void setLevelCouplings() {}
void setBoundaryConditions();
auto& weakPatches() {
return weakPatches_;
}
private:
static constexpr Dune::FieldVector<double, MY_DIM> zenith_() {
#if MY_DIM == 2
return {0, 1};
#elif MY_DIM == 3
return {0, 1, 0};
#endif
}
};
#endif
dune/tectonic/factories/twoblocksfactory_tmpl.cc 0 → 100644
View file @ 89161286
#ifndef MY_DIM
#error MY_DIM unset
#endif
#include "../explicitgrid.hh"
#include "../explicitvectors.hh"
template class TwoBlocksFactory<Grid, Vector>;
dune/tectonic/globalratestatefriction.hh deleted 100644 → 0
View file @ 8fe950bc
#ifndef DUNE_TECTONIC_GLOBALRATESTATEFRICTION_HH
#define DUNE_TECTONIC_GLOBALRATESTATEFRICTION_HH
#include <vector>
#include <dune/common/bitsetvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/tectonic/geocoordinate.hh>
#include <dune/tectonic/globalfrictiondata.hh>
#include <dune/tectonic/globalfriction.hh>
#include <dune/tectonic/index-in-sorted-range.hh>
template <class Matrix, class Vector, class ScalarFriction, class GridView>
class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
public:
using GlobalFriction<Matrix, Vector>::block_size;
using typename GlobalFriction<Matrix, Vector>::LocalNonlinearity;
private:
using typename GlobalFriction<Matrix, Vector>::ScalarVector;
public:
GlobalRateStateFriction(BoundaryPatch<GridView> const &frictionalBoundary,
GlobalFrictionData<block_size> const &frictionInfo,
ScalarVector const &weights,
ScalarVector const &weightedNormalStress)
: restrictions(), localToGlobal(), zeroFriction() {
auto const gridView = frictionalBoundary.gridView();
Dune::MultipleCodimMultipleGeomTypeMapper<
GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView);
for (auto it = gridView.template begin<block_size>();
it != gridView.template end<block_size>(); ++it) {
auto const i = vertexMapper.index(*it);
if (not frictionalBoundary.containsVertex(i))
continue;
localToGlobal.emplace_back(i);
restrictions.emplace_back(weights[i], weightedNormalStress[i],
frictionInfo(geoToPoint(it->geometry())));
}
assert(restrictions.size() == frictionalBoundary.numVertices());
assert(localToGlobal.size() == frictionalBoundary.numVertices());
}
void updateAlpha(ScalarVector const &alpha) override {
for (size_t j = 0; j < restrictions.size(); ++j)
restrictions[j].updateAlpha(alpha[localToGlobal[j]]);
}
/*
Return a restriction of the outer function to the i'th node.
*/
LocalNonlinearity const &restriction(size_t i) const override {
auto const index = indexInSortedRange(localToGlobal, i);
if (index == localToGlobal.size())
return zeroFriction;
return restrictions[index];
}
private:
std::vector<WrappedScalarFriction<block_size, ScalarFriction>> restrictions;
std::vector<size_t> localToGlobal;
WrappedScalarFriction<block_size, ZeroFunction> const zeroFriction;
};
#endif
src/gridselector.hh → dune/tectonic/gridselector.hh
View file @ 89161286
File moved
dune/tectonic/io/CMakeLists.txt 0 → 100644
View file @ 89161286
add_subdirectory("hdf5")
add_custom_target(tectonic_dune_io SOURCES
hdf5-bodywriter.hh
hdf5-levelwriter.hh
uniform-grid-writer.cc
vtk.hh
vtk.cc
)
#install headers
install(FILES
hdf5-bodywriter.hh
hdf5-levelwriter.hh
vtk.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)
dune/tectonic/io/hdf5-bodywriter.hh 0 → 100644
View file @ 89161286
#ifndef SRC_HDF5_BODYWRITER_HH
#define SRC_HDF5_BODYWRITER_HH
#include <dune/fufem/functions/basisgridfunction.hh>
#include <dune/fufem/geometry/convexpolyhedron.hh>
#include <dune/fufem/hdf5/file.hh>
#include "hdf5/frictionalboundary-writer.hh"
#include "hdf5/patchinfo-writer.hh"
template <class ProgramState, class VertexBasis, class GridView>
class HDF5BodyWriter {
private:
using PatchInfoWriter = PatchInfoWriter<ProgramState, VertexBasis, GridView>;
using FrictionalBoundaryWriter = FrictionalBoundaryWriter<ProgramState, GridView>;
public:
using VertexCoordinates = typename ProgramState::Vector;
using FrictionPatches = std::vector<const BoundaryPatch<GridView>* >;
using WeakPatches = std::vector<const ConvexPolyhedron<LocalVector>* >;
using LocalVector = typename VertexCoordinates::block_type;
//friend class HDF5LevelWriter<ProgramState, VertexBasis, GridView>;
HDF5BodyWriter(
HDF5::Grouplike &file,
const VertexCoordinates& vertexCoordinates,
const VertexBasis& vertexBasis,
const FrictionPatches& frictionPatches,
const WeakPatches& weakPatches) :
patchCount_(frictionPatches.size()),
file_(file),
#if MY_DIM == 3
patchInfoWriters_(patchCount_),
#endif
frictionBoundaryWriters_(patchCount_) {
assert(patchCount_ == weakPatches.size());
for (size_t i=0; i<patchCount_; i++) {
#if MY_DIM == 3
patchInfoWriters_[i] = std::make_unique<PatchInfoWriter>(file_, vertexBasis, *frictionPatches[i], *weakPatches[i], i);
#endif
frictionBoundaryWriters_[i] = std::make_unique<FrictionalBoundaryWriter>(file_, vertexCoordinates, *frictionPatches[i], i);
}
}
template <class Friction>
void reportSolution(ProgramState const &programState,
// for the friction coefficient
Friction& friction) {
for (size_t i=0; i<patchCount_; i++) {
#if MY_DIM == 3
patchInfoWriters_[i]->write(programState);
#endif
frictionBoundaryWriters_[i]->write(programState, friction);
}
}
private:
const size_t patchCount_;
HDF5::Grouplike &file_;
#if MY_DIM == 3
std::vector<std::unique_ptr<PatchInfoWriter>> patchInfoWriters_;
#endif
std::vector<std::unique_ptr<FrictionalBoundaryWriter>> frictionBoundaryWriters_;
};
#endif
dune/tectonic/io/hdf5-levelwriter.hh 0 → 100644
View file @ 89161286
#ifndef SRC_HDF5_LEVELWRITER_HH
#define SRC_HDF5_LEVELWRITER_HH
#include <dune/fufem/hdf5/file.hh>
#include "hdf5/iteration-writer.hh"
#include "hdf5/time-writer.hh"
#include "hdf5-bodywriter.hh"
template <class ProgramState, class VertexBasis, class GridView>
class HDF5LevelWriter {
public:
using HDF5BodyWriter = HDF5BodyWriter<ProgramState, VertexBasis, GridView>;
using VertexCoordinates = std::vector<typename HDF5BodyWriter::VertexCoordinates>;
using VertexBases = std::vector<const VertexBasis*>;
using FrictionPatches = std::vector<typename HDF5BodyWriter::FrictionPatches>;
using WeakPatches = std::vector<typename HDF5BodyWriter::WeakPatches>;
//friend class HDF5NetworkWriter<ProgramState, VertexBasis, GridView>;
HDF5LevelWriter(HDF5::Grouplike &file,
const VertexCoordinates& vertexCoordinates,
const VertexBases& vertexBases,
const FrictionPatches& frictionBoundaries,
const WeakPatches& weakPatches)
: bodyCount_(vertexCoordinates.size()),
file_(file),
iterationWriter_(file_),
timeWriter_(file_),
bodyWriters_(bodyCount_) {
for (size_t i=0; i<bodyCount_; i++) {
bodyWriters_[i] = std::make_unique<HDF5BodyWriter>(file_, vertexCoordinates[i], *vertexBases[i], frictionBoundaries[i], weakPatches[i]);
}
}
template <class Friction>
void reportSolution(
const ProgramState& programState,
// for the friction coefficient
Friction& friction) {
timeWriter_.write(programState);
for (size_t i=0; i<bodyCount_; i++) {
bodyWriters_[i]->reportSolution(programState, friction); //TODO
}
}
void reportIterations(
const ProgramState& programState,
const IterationRegister& iterationCount) {
iterationWriter_.write(programState.timeStep, iterationCount);
}
private:
const size_t bodyCount_;
HDF5::Grouplike &file_;
IterationWriter iterationWriter_;
TimeWriter<ProgramState> timeWriter_;
std::vector<std::unique_ptr<HDF5BodyWriter>> bodyWriters_;
};
#endif
dune/tectonic/io/hdf5/CMakeLists.txt 0 → 100644
View file @ 89161286
add_custom_target(tectonic_dune_io_hdf5 SOURCES
frictionalboundary-writer.hh
frictionalboundary-writer.cc
iteration-writer.hh
iteration-writer.cc
patchinfo-writer.hh
patchinfo-writer.cc
restart-io.hh
restart-io.cc
restrict.hh
surface-writer.hh
surface-writer.cc
time-writer.hh
time-writer.cc
)
#install headers
install(FILES
frictionalboundary-writer.hh
iteration-writer.hh
patchinfo-writer.hh
restart-io.hh
restrict.hh
surface-writer.hh
time-writer.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)
src/hdf5/frictionalboundary-writer.cc → dune/tectonic/io/hdf5/frictionalboundary-writer.cc
View file @ 89161286
...@@ -10,8 +10,9 @@ ...@@ -10,8 +10,9 @@
template <class ProgramState, class GridView> template <class ProgramState, class GridView>
FrictionalBoundaryWriter<ProgramState, GridView>::FrictionalBoundaryWriter( FrictionalBoundaryWriter<ProgramState, GridView>::FrictionalBoundaryWriter(
HDF5::Grouplike &file, Vector const &vertexCoordinates, HDF5::Grouplike &file, Vector const &vertexCoordinates,
Patch const &frictionalBoundary) Patch const &frictionalBoundary, size_t const id)
: group_(file, "frictionalBoundary"), : id_(id),
group_(file, "frictionalBoundary"+std::to_string(id_)),
frictionalBoundary_(frictionalBoundary), frictionalBoundary_(frictionalBoundary),
frictionalBoundaryDisplacementWriter_(group_, "displacement", frictionalBoundaryDisplacementWriter_(group_, "displacement",
frictionalBoundary.numVertices(), frictionalBoundary.numVertices(),
...@@ -35,27 +36,29 @@ template <class ProgramState, class GridView> ...@@ -35,27 +36,29 @@ template <class ProgramState, class GridView>
template <class Friction> template <class Friction>
void FrictionalBoundaryWriter<ProgramState, GridView>::write( void FrictionalBoundaryWriter<ProgramState, GridView>::write(
ProgramState const &programState, Friction &friction) { ProgramState const &programState, Friction &friction) {
auto const frictionalBoundaryDisplacements = auto const frictionalBoundaryDisplacements =
restrictToSurface(programState.u, frictionalBoundary_); restrictToSurface(programState.u[id_], frictionalBoundary_);
addEntry(frictionalBoundaryDisplacementWriter_, programState.timeStep, addEntry(frictionalBoundaryDisplacementWriter_, programState.timeStep,
frictionalBoundaryDisplacements); frictionalBoundaryDisplacements);
auto const frictionalBoundaryVelocities = auto const frictionalBoundaryVelocities =
restrictToSurface(programState.v, frictionalBoundary_); restrictToSurface(programState.v[id_], frictionalBoundary_);
addEntry(frictionalBoundaryVelocityWriter_, programState.timeStep, addEntry(frictionalBoundaryVelocityWriter_, programState.timeStep,
frictionalBoundaryVelocities); frictionalBoundaryVelocities);
auto const frictionalBoundaryStates = auto const frictionalBoundaryStates =
restrictToSurface(programState.alpha, frictionalBoundary_); restrictToSurface(programState.alpha[id_], frictionalBoundary_);
addEntry(frictionalBoundaryStateWriter_, programState.timeStep, addEntry(frictionalBoundaryStateWriter_, programState.timeStep,
frictionalBoundaryStates); frictionalBoundaryStates);
friction.updateAlpha(programState.alpha); //TODO: needs to transform programState.v to relative velocities with nBodyAssembler
auto const c = friction.coefficientOfFriction(programState.v); /*friction.updateAlpha(programState.alpha);
auto const c = friction.coefficientOfFriction(programState.v[id_]);
auto const frictionalBoundaryCoefficient = auto const frictionalBoundaryCoefficient =
restrictToSurface(c, frictionalBoundary_); restrictToSurface(c, frictionalBoundary_);
addEntry(frictionalBoundaryCoefficientWriter_, programState.timeStep, addEntry(frictionalBoundaryCoefficientWriter_, programState.timeStep,
frictionalBoundaryCoefficient); frictionalBoundaryCoefficient);*/
} }
#include "frictionalboundary-writer_tmpl.cc" #include "frictionalboundary-writer_tmpl.cc"
src/hdf5/frictionalboundary-writer.hh → dune/tectonic/io/hdf5/frictionalboundary-writer.hh
View file @ 89161286
...@@ -11,13 +11,15 @@ template <class ProgramState, class GridView> class FrictionalBoundaryWriter { ...@@ -11,13 +11,15 @@ template <class ProgramState, class GridView> class FrictionalBoundaryWriter {
using Patch = BoundaryPatch<GridView>; using Patch = BoundaryPatch<GridView>;
public: public:
FrictionalBoundaryWriter(HDF5::Grouplike &file, Vector const &vertexCoordinates, FrictionalBoundaryWriter(HDF5::Grouplike& file, Vector const &vertexCoordinates,
Patch const &frictionalBoundary); Patch const &frictionalBoundary, size_t const id);
template <class Friction> template <class Friction>
void write(ProgramState const &programState, Friction &friction); void write(ProgramState const &programState, Friction &friction);
private: private:
size_t const id_;
HDF5::Group group_; HDF5::Group group_;
Patch const &frictionalBoundary_; Patch const &frictionalBoundary_;
......
src/hdf5/frictionalboundary-writer_tmpl.cc → dune/tectonic/io/hdf5/frictionalboundary-writer_tmpl.cc
View file @ 89161286
#include "../explicitvectors.hh" #include "../../explicitvectors.hh"
#include "../explicitgrid.hh" #include "../../explicitgrid.hh"
#include "../program_state.hh" #include "../../data-structures/program_state.hh"
using MyProgramState = ProgramState<Vector, ScalarVector>; using MyProgramState = ProgramState<Vector, ScalarVector>;
using MyFriction = GlobalFriction<Matrix, Vector>; using MyFriction = GlobalFriction<Matrix, Vector>;
template class FrictionalBoundaryWriter<MyProgramState, GridView>; template class FrictionalBoundaryWriter<MyProgramState, DeformedGrid::LeafGridView>;
template void FrictionalBoundaryWriter<MyProgramState, GridView>::write( template void FrictionalBoundaryWriter<MyProgramState, DeformedGrid::LeafGridView>::write(
MyProgramState const &programState, MyFriction &friction); MyProgramState const &programState, MyFriction &friction);