From 0a89a40f89eee607b047c68873964e01d2f6f8aa Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@mi.fu-berlin.de>
Date: Thu, 11 Apr 2019 18:10:59 +0200
Subject: [PATCH] .
---
src/CMakeLists.txt | 1 +
src/nodalweights.cc | 75 +++++++
src/nodalweights.hh | 37 ++++
src/tests/CMakeLists.txt | 2 +
src/tests/common.hh | 62 ++++++
src/tests/gridgluefrictiontest.cc | 327 ++++++++++++++++++++++++++++++
src/tests/nodalweightstest.cc | 170 ++++++++++++++++
src/utils/debugutils.hh | 98 ++++++++-
8 files changed, 762 insertions(+), 10 deletions(-)
create mode 100644 src/nodalweights.cc
create mode 100644 src/nodalweights.hh
create mode 100644 src/tests/common.hh
create mode 100644 src/tests/gridgluefrictiontest.cc
create mode 100644 src/tests/nodalweightstest.cc
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 8a218775..32c35011 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -25,6 +25,7 @@ set(SW_SOURCE_FILES
set(MSW_SOURCE_FILES
assemblers.cc
+ nodalweights.cc
data-structures/body.cc
#data-structures/contactnetwork.cc
data-structures/enumparser.cc
diff --git a/src/nodalweights.cc b/src/nodalweights.cc
new file mode 100644
index 00000000..910585f1
--- /dev/null
+++ b/src/nodalweights.cc
@@ -0,0 +1,75 @@
+#include "utils/almostequal.hh"
+
+#include "nodalweights.hh"
+
+template <class Basis0, class Basis1>
+NodalWeights<Basis0, Basis1>::NodalWeights(const Basis0& basis0, const Basis1& basis1)
+ : basis0_(basis0),
+ basis1_(basis1)
+ {}
+
+template <class Basis0, class Basis1>
+template <class Basis, class Element, class GlobalCoords>
+auto NodalWeights<Basis0, Basis1>::basisDof(const Basis& basis, const Element& elem, const GlobalCoords& vertex) const {
+ int dof = -1;
+
+ const typename Basis::LocalFiniteElement& lFE = basis.getLocalFiniteElement(elem);
+ const auto& geometry = elem.geometry();
+ const auto& localVertex = geometry.local(vertex);
+
+ const size_t localBasisSize = lFE.localBasis().size();
+ std::vector<Dune::FieldVector<double, 1>, std::allocator<Dune::FieldVector<double, 1> > > localBasisRep(localBasisSize);
+ lFE.localBasis().evaluateFunction(localVertex, localBasisRep);
+
+ for(size_t i=0; i<localBasisSize; i++) {
+ if (almost_equal(localBasisRep[i][0], 1.0, 2)) {
+ dof = basis.index(elem, i);
+ break;
+ }
+ }
+
+ return dof;
+}
+
+template <class Basis0, class Basis1>
+template <class GridGlue, class ScalarVector>
+void NodalWeights<Basis0, Basis1>::assemble(const GridGlue& glue, ScalarVector& weights0, ScalarVector& weights1, bool initializeVector) const {
+ using ctype = typename ScalarVector::field_type;
+
+ if (initializeVector==true) {
+ weights0.resize(basis0_.size());
+ weights1.resize(basis1_.size());
+ }
+
+ // loop over all intersections
+ for (const auto& rIs : intersections(glue)) {
+ const auto& inside = rIs.inside();
+ const auto& outside = rIs.outside();
+
+ /*if (!nmBoundary.contains(inside, rIs.indexInInside())) {
+ std::cout << "it happened" << std::endl;
+ continue;
+ }*/
+
+ const auto& geometry = rIs.geometry();
+ const ctype val = 1.0/(geometry.mydimension + 1)*geometry.volume();
+
+ std::cout << geometry.mydimension << " " << geometry.volume() << " " << val << std::endl;
+
+ for (int i=0; i<geometry.corners(); i++) {
+ const auto& vertex = geometry.corner(i);
+
+ const int inIdx = basisDof(basis0_, inside, vertex);
+ if (inIdx >= 0) {
+ weights0[inIdx] += val;
+ }
+
+ const int outIdx = basisDof(basis1_, outside, vertex);
+ if (outIdx >= 0) {
+ weights1[outIdx] += val;
+ }
+
+ std::cout << inIdx << " " << outIdx << std::endl;
+ }
+ }
+}
diff --git a/src/nodalweights.hh b/src/nodalweights.hh
new file mode 100644
index 00000000..bafd4d96
--- /dev/null
+++ b/src/nodalweights.hh
@@ -0,0 +1,37 @@
+#ifndef SRC_NODALWEIGHTS_HH
+#define SRC_NODALWEIGHTS_HH
+
+
+/**
+ * let merged contact boundary Gamma be given by GridGlue object;
+ * computes
+ * int_Gamma lambda_i dx,
+ * where lambda_i is nodal basis of merged contact boundary whose
+ * dofs are given by nonmortar and mortar dofs;
+ *
+ *
+ * NOTE: works only for P1 nodal bases
+ **/
+
+template <class Basis0, class Basis1>
+class NodalWeights {
+private:
+ enum {dim = Basis0::GridView::Grid::dimension};
+
+ template <class Basis, class Element, class GlobalCoords>
+ auto basisDof(const Basis& basis, const Element& elem, const GlobalCoords& vertex) const;
+
+public:
+ NodalWeights(const Basis0& basis0, const Basis1& basis1);
+
+ template <class GridGlue, class ScalarVector>
+ void assemble(const GridGlue& glue, ScalarVector& weights0, ScalarVector& weights1, bool initializeVector = true) const;
+
+private:
+ const Basis0& basis0_;
+ const Basis1& basis1_;
+};
+
+#include "nodalweights.cc"
+
+#endif
diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt
index 4109f1b9..8ca00f64 100644
--- a/src/tests/CMakeLists.txt
+++ b/src/tests/CMakeLists.txt
@@ -1 +1,3 @@
dune_add_test(SOURCES globalfrictioncontainertest.cc)
+dune_add_test(SOURCES gridgluefrictiontest.cc)
+dune_add_test(SOURCES nodalweightstest.cc)
diff --git a/src/tests/common.hh b/src/tests/common.hh
new file mode 100644
index 00000000..7e679834
--- /dev/null
+++ b/src/tests/common.hh
@@ -0,0 +1,62 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_TECTONIC_SRC_TESTS_COMMON_HH
+#define DUNE_TECTONIC_SRC_TESTS_COMMON_HH
+
+#include <dune/grid/common/gridfactory.hh>
+#include <dune/grid/utility/structuredgridfactory.hh>
+
+#include <dune/fufem/referenceelementhelper.hh>
+
+// utility functions
+bool isClose(double a, double b) {
+ return std::abs(a - b) < 1e-14;
+}
+
+template <class Vector>
+bool xyCollinear(Vector const &a, Vector const &b, Vector const &c) {
+ return isClose((b[0] - a[0]) * (c[1] - a[1]), (b[1] - a[1]) * (c[0] - a[0]));
+}
+
+template <class Vector>
+bool xyBoxed(Vector const &v1, Vector const &v2, Vector const &x) {
+ auto const minmax0 = std::minmax(v1[0], v2[0]);
+ auto const minmax1 = std::minmax(v1[1], v2[1]);
+
+ if (minmax0.first - 1e-14 > x[0] or
+ x[0] > minmax0.second + 1e-14)
+ return false;
+ if (minmax1.first - 1e-14 > x[1] or
+ x[1] > minmax1.second + 1e-14)
+ return false;
+
+ return true;
+}
+
+template <class Vector>
+bool xyBetween(Vector const &v1, Vector const &v2, Vector const &x) {
+ return xyCollinear(v1, v2, x) && xyBoxed(v1, v2, x);
+}
+
+// build cube grid given by lowerLeft and upperRight point in global coordinates
+template <class GlobalCoords, class GridType>
+void buildGrid(const GlobalCoords& lowerLeft, const GlobalCoords& upperRight, const int n, std::shared_ptr<GridType>& grid, bool simplexGrid = true) {
+ std::array<unsigned int, GridType::dimension> elements;
+ std::fill(elements.begin(), elements.end(), n);
+
+ Dune::GridFactory<GridType> factory;
+ if (!simplexGrid) {
+ Dune::StructuredGridFactory<GridType>::createCubeGrid(factory, lowerLeft, upperRight, elements);
+ grid = std::move(factory.createGrid());
+ } else {
+ Dune::StructuredGridFactory<GridType>::createSimplexGrid(factory, lowerLeft, upperRight, elements);
+ grid = std::move(factory.createGrid());
+ }
+}
+
+template <class Intersection>
+bool containsInsideSubentity(const Intersection& nIt, int subEntity, int codim) {
+ return ReferenceElementHelper<double, Intersection::GridView::dim>::subEntityContainsSubEntity(nIt.inside().type(), nIt.indexInInside(), 1, subEntity, codim);
+}
+
+#endif
diff --git a/src/tests/gridgluefrictiontest.cc b/src/tests/gridgluefrictiontest.cc
new file mode 100644
index 00000000..e66fa3dd
--- /dev/null
+++ b/src/tests/gridgluefrictiontest.cc
@@ -0,0 +1,327 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <exception>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/stdstreams.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/function.hh>
+
+//#include <dune/fufem/assemblers/assembler.hh>
+#include <dune/fufem/functiontools/basisinterpolator.hh>
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+
+#include <dune/fufem/functions/basisgridfunction.hh>
+
+#include <dune/grid/uggrid.hh>
+
+#include <dune/grid-glue/adapter/gridgluevtkwriter.hh>
+#include <dune/grid-glue/merging/merger.hh>
+
+#include <dune/contact/projections/normalprojection.hh>
+#include <dune/contact/common/couplingpair.hh>
+#include <dune/contact/assemblers/dualmortarcoupling.hh>
+#include <dune/contact/assemblers/nbodyassembler.hh>
+
+#include "../utils/debugutils.hh"
+
+#include "common.hh"
+
+const int dim = 2;
+const int n = 5;
+const bool simplexGrid = true;
+
+const std::string path = "";
+const std::string outputFile = "gridgluefrictiontest.log";
+
+#if HAVE_UG
+ using GridType = typename Dune::UGGrid<dim>;
+#else
+#error No UG!
+#endif
+ using LevelView = typename GridType::LevelGridView;
+ using LeafView = typename GridType::LeafGridView;
+ using LevelBoundaryPatch = BoundaryPatch<LevelView>;
+ using LeafBoundaryPatch = BoundaryPatch<LeafView>;
+
+ using c_type = double;
+ using GlobalCoords = Dune::FieldVector<c_type, dim>;
+ using Vector = Dune::FieldVector<c_type, 1>;
+ using BlockVector = Dune::BlockVector<Vector>;
+
+ using CouplingPair = Dune::Contact::CouplingPair<GridType, GridType, c_type>;
+ using CouplingType = Dune::Contact::DualMortarCoupling<c_type, GridType>;
+ using NBodyAssembler = Dune::Contact::NBodyAssembler<GridType, BlockVector>;
+
+
+
+// boundary function along y=1.0 line
+template <class VectorType>
+class BoundaryFunction0 : public Dune::VirtualFunction<VectorType, c_type> {
+public:
+ void evaluate(const VectorType& x, c_type& res) const override {
+ if (isClose(x[1], 1.0)) {
+ res = std::sin(x[0]);
+ } else {
+ res = 0.0;
+ }
+ }
+};
+
+// boundary function along y=1.0 line
+template <class VectorType>
+class BoundaryFunction1 : public Dune::VirtualFunction<VectorType, c_type> {
+public:
+ void evaluate(const VectorType& x, c_type& res) const override{
+ if (isClose(x[1], 1.0)) {
+ res = std::cos(x[0]);
+ } else {
+ res = 0.0;
+ }
+ }
+};
+
+
+int main(int argc, char *argv[]) { try {
+ Dune::MPIHelper::instance(argc, argv);
+
+ std::ofstream out(path + outputFile);
+ std::streambuf *coutbuf = std::cout.rdbuf(); //save old buffer
+ std::cout.rdbuf(out.rdbuf()); //redirect std::cout to outputFile
+
+ std::cout << "------------------------------" << std::endl;
+ std::cout << "-- Grid-Glue Friction Test: --" << std::endl;
+ std::cout << "------------------------------" << std::endl << std::endl;
+
+ // building grids
+ std::vector<std::shared_ptr<GridType>> grids(2);
+
+ GlobalCoords lowerLeft0({0, 0});
+ GlobalCoords upperRight0({2, 1});
+ buildGrid(lowerLeft0, upperRight0, n, grids[0]);
+
+ GlobalCoords lowerLeft1({0.25, 1});
+ GlobalCoords upperRight1({1.75, 2});
+ buildGrid(lowerLeft1, upperRight1, n, grids[1]);
+
+ // writing grids
+ for (size_t i=0; i<grids.size(); i++) {
+ const auto& levelView = grids[i]->levelGridView(0);
+ writeToVTK(levelView, path, "body_" + std::to_string(i) + "_level0");
+ }
+
+ // compute coupling boundaries
+ LevelView gridView0 = grids[0]->levelGridView(0);
+ LevelView gridView1 = grids[1]->levelGridView(0);
+ LevelBoundaryPatch upper(gridView0);
+ LevelBoundaryPatch lower(gridView1);
+
+ lower.insertFacesByProperty([&](typename LevelView::Intersection const &in) {
+ return xyBetween(lowerLeft1, {upperRight1[0], lowerLeft1[1]}, in.geometry().center());
+ });
+
+ upper.insertFacesByProperty([&](typename LevelView::Intersection const &in) {
+ return xyBetween({lowerLeft0[0], upperRight0[1]}, upperRight0, in.geometry().center());
+ });
+
+ // set contact coupling
+ Dune::Contact::NormalProjection<LeafBoundaryPatch> contactProjection;
+
+ CouplingPair coupling;
+ coupling.set(0, 1, upper, lower, 0.1, CouplingPair::CouplingType::STICK_SLIP, contactProjection, nullptr);
+
+ double coveredArea_ = 0.8;
+ CouplingType contactCoupling;
+ contactCoupling.setGrids(*grids[0], *grids[1]);
+ contactCoupling.setupContactPatch(*coupling.patch0(),*coupling.patch1());
+ contactCoupling.gridGlueBackend_ = coupling.backend();
+ contactCoupling.setCoveredArea(coveredArea_);
+ contactCoupling.setup();
+
+ // set nBodyAssembler
+ /* std::vector<const GridType*> grids_ptr(grids.size());
+ for (size_t i=0; i<grids_ptr.size(); i++) {
+ grids_ptr[i] = grids[i].get();
+ }
+
+ NBodyAssembler nBodyAssembler(grids.size(), 1);
+ nBodyAssembler.setGrids(grids_ptr);
+
+ nBodyAssembler.setCoupling(coupling, 0);
+
+ nBodyAssembler.assembleTransferOperator();
+ nBodyAssembler.assembleObstacle(); */
+
+ // define basis
+ using Basis = P1NodalBasis<LevelView, double>;
+ Basis basis0(grids[0]->levelGridView(0));
+ Basis basis1(grids[1]->levelGridView(0));
+
+ std::cout << "--------------" << std::endl;
+ std::cout << "--- Basis0 ---" << std::endl;
+ std::cout << "--------------" << std::endl;
+ printBasisDofLocation(basis0);
+
+ std::cout << "--------------" << std::endl;
+ std::cout << "--- Basis1 ---" << std::endl;
+ std::cout << "--------------" << std::endl;
+ printBasisDofLocation(basis1);
+
+ // set grid functions on coupling boundary
+ BlockVector f0, f1;
+
+ BasisInterpolator<Basis, BlockVector, BoundaryFunction0<GlobalCoords>> interpolator0;
+ interpolator0.interpolate(basis0, f0, BoundaryFunction0<GlobalCoords>());
+
+ BasisInterpolator<Basis, BlockVector, BoundaryFunction1<GlobalCoords>> interpolator1;
+ interpolator1.interpolate(basis1, f1, BoundaryFunction1<GlobalCoords>());
+
+ printBoundary(upper, basis0, f0, "f0: ");
+ printBoundary(lower, basis1, f1, "f1: ");
+ writeToVTK(basis0, f0, path, "body_0_level0");
+ writeToVTK(basis1, f1, path, "body_1_level0");
+
+ // merged gridGlue coupling boundary
+ auto& gridGlue = *contactCoupling.getGlue();
+
+ // make basis grid functions
+ auto&& gridFunction0 = Functions::makeFunction(basis0, f0);
+ auto&& gridFunction1 = Functions::makeFunction(basis1, f1);
+
+ for (const auto& rIs : intersections(gridGlue)) {
+ const auto& inside = rIs.inside();
+ const auto& outside = rIs.outside();
+
+ const auto& rGeometry = rIs.geometry();
+ const auto& outGeometry = outside.geometry();
+
+ for (size_t i=0; i<rGeometry.corners(); i++) {
+ typename BasisGridFunction<Basis, BlockVector>::RangeType y;
+ gridFunction1.evaluateLocal(outside, outGeometry.local(rGeometry.corner(i)), y);
+ print(rGeometry.corner(i), "corner " + std::to_string(i));
+ std::cout << "f1(corner) = " << y << std::endl;
+
+ std::cout << std::endl;
+ }
+ std::cout << "---------"<< std::endl;
+ // auto basis.
+ // std::cout << containsInsideSubentity(rIs, 1, 1) << std::endl;
+ }
+
+ //Dune::GridGlue::GridGlueVtkWriter::writeIntersections<std::decltype(gridGlue), 0>(gridGlue, "");//, path + "gridglueintersections0.vtk"); //<std::decltype(gridGlue), 0>
+
+ /*LevelBoundaryPatch boundaryWithMapping(grids[0]->levelGridView(0));
+
+ const auto& indexSet0 = gridView0.indexSet();
+ const auto& indexSet1 = gridView1.indexSet();
+
+ // Get all fine grid boundary segments that are totally covered by the grid-glue segments
+ typedef std::pair<int,int> Pair;
+ std::map<Pair, c_type> coveredArea, fullArea;
+
+ // initialize with area of boundary faces
+ for (const auto& bIt : lower) {
+ const Pair p(indexSet0.index(bIt.inside()),bIt.indexInInside());
+ fullArea[p] = bIt.geometry().volume();
+ coveredArea[p] = 0;
+ }
+
+ // sum up the remote intersection areas to find out which are totally covered
+ for (const auto& rIs : intersections(gridGlue))
+ coveredArea[Pair(indexSet0.index(rIs.inside()),rIs.indexInInside())] += rIs.geometry().volume();
+
+ // add all fine grid faces that are totally covered by the contact mapping
+ for (const auto& bIt : lower) {
+ const auto& inside = bIt.inside();
+ if(coveredArea[Pair(indexSet0.index(inside),bIt.indexInInside())]/
+ fullArea[Pair(indexSet0.index(inside),bIt.indexInInside())] >= coveredArea_)
+ boundaryWithMapping.addFace(inside, bIt.indexInInside());
+ }*/
+
+ //writeBoundary(boundaryWithMapping, path + "relevantNonmortar");
+
+
+ /** \todo replace by all fine grid segments which are totally covered by the RemoteIntersections. */
+ //for (const auto& rIs : intersections(glue))
+ // boundaryWithMapping.addFace(rIs.inside(),rIs.indexInInside());
+
+ /* printf("contact mapping could be built for %d of %d boundary segments.\n",
+ boundaryWithMapping.numFaces(), lower.numFaces());
+
+ LevelBoundaryPatch nonmortarBoundary = boundaryWithMapping;
+ LevelBoundaryPatch mortarBoundary;
+ mortarBoundary.setup(gridView1);
+ for (const auto& rIs : intersections(gridGlue))
+ if (nonmortarBoundary.contains(rIs.inside(),rIs.indexInInside()))
+ mortarBoundary.addFace(rIs.outside(),rIs.indexInOutside());*/
+
+
+
+
+
+
+ /** \todo Also restrict mortar indices and don't use the whole grid level. */
+ /* Dune::MatrixIndexSet mortarIndices(nonmortarBoundary_.numVertices(), grid1_->size(dim));
+
+ // Create mapping from the global set of block dofs to the ones on the contact boundary
+ std::vector<int> globalToLocal;
+ nonmortarBoundary_.makeGlobalToLocal(globalToLocal);
+
+ // loop over all intersections
+ for (const auto& rIs : intersections(glue)) {
+
+ if (!nonmortarBoundary_.contains(rIs.inside(),rIs.indexInInside()))
+ continue;
+
+ const auto& inside = rIs.inside();
+ const auto& outside = rIs.outside();
+
+ const auto& domainRefElement = Dune::ReferenceElements<ctype, dim>::general(inside.type());
+ const auto& targetRefElement = Dune::ReferenceElements<ctype, dim>::general(outside.type());
+
+ int nDomainVertices = domainRefElement.size(dim);
+ int nTargetVertices = targetRefElement.size(dim);
+
+ for (int j=0; j<nDomainVertices; j++) {
+
+ int localDomainIdx = globalToLocal[indexSet0.subIndex(inside,j,dim)];
+
+ // if the vertex is not contained in the restricted contact boundary then dismiss it
+ if (localDomainIdx == -1)
+ continue;
+
+ for (int k=0; k<nTargetVertices; k++) {
+ int globalTargetIdx = indexSet1.subIndex(outside,k,dim);
+ if (!mortarBoundary_.containsVertex(globalTargetIdx))
+ continue;
+
+ mortarIndices.add(localDomainIdx, globalTargetIdx);
+ }
+ }
+ }
+
+ mortarIndices.exportIdx(mortarLagrangeMatrix_);
+
+ // Clear it
+ mortarLagrangeMatrix_ = 0; */
+
+ bool passed = true;
+
+ std::cout << "Overall, the test " << (passed ? "was successful!" : "failed!") << std::endl;
+
+ std::cout.rdbuf(coutbuf); //reset to standard output again
+ return passed ? 0 : 1;
+
+} catch (Dune::Exception &e) {
+ Dune::derr << "Dune reported error: " << e << std::endl;
+} catch (std::exception &e) {
+ std::cerr << "Standard exception: " << e.what() << std::endl;
+} // end try
+} // end main
diff --git a/src/tests/nodalweightstest.cc b/src/tests/nodalweightstest.cc
new file mode 100644
index 00000000..78b3c7c2
--- /dev/null
+++ b/src/tests/nodalweightstest.cc
@@ -0,0 +1,170 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <iostream>
+#include <exception>
+#include <fstream>
+
+#include <dune/grid/uggrid.hh>
+
+#include <dune/fufem/assemblers/assembler.hh>
+#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
+#include <dune/fufem/functions/constantfunction.hh>
+#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+
+#include <dune/contact/projections/normalprojection.hh>
+#include <dune/contact/common/couplingpair.hh>
+#include <dune/contact/assemblers/dualmortarcoupling.hh>
+
+#include "../utils/debugutils.hh"
+#include "common.hh"
+
+#include "../nodalweights.hh"
+
+const int dim = 2;
+const int n = 5;
+const bool simplexGrid = true;
+
+const std::string path = "";
+const std::string outputFile = "nodalweightstest.log";
+
+#if HAVE_UG
+ using GridType = typename Dune::UGGrid<dim>;
+#else
+#error No UG!
+#endif
+ using LevelView = typename GridType::LevelGridView;
+ using LeafView = typename GridType::LeafGridView;
+ using LevelBoundaryPatch = BoundaryPatch<LevelView>;
+ using LeafBoundaryPatch = BoundaryPatch<LeafView>;
+
+ using c_type = double;
+ using GlobalCoords = Dune::FieldVector<c_type, dim>;
+ using BlockVector = Dune::BlockVector<GlobalCoords>;
+
+ using ScalarVector = Dune::FieldVector<c_type, 1>;
+ using ScalarBlockVector = Dune::BlockVector<ScalarVector>;
+
+ using CouplingPair = Dune::Contact::CouplingPair<GridType, GridType, c_type>;
+ using CouplingType = Dune::Contact::DualMortarCoupling<c_type, GridType>;
+
+int main(int argc, char *argv[]) { try {
+ Dune::MPIHelper::instance(argc, argv);
+
+ std::ofstream out(path + outputFile);
+ std::streambuf *coutbuf = std::cout.rdbuf(); //save old buffer
+ std::cout.rdbuf(out.rdbuf()); //redirect std::cout to outputFile
+
+ std::cout << "------------------------" << std::endl;
+ std::cout << "-- NodalWeights Test: --" << std::endl;
+ std::cout << "------------------------" << std::endl << std::endl;
+
+ // building grids
+ std::vector<std::shared_ptr<GridType>> grids(2);
+
+ GlobalCoords lowerLeft0({0, 0});
+ GlobalCoords upperRight0({2, 1});
+ buildGrid(lowerLeft0, upperRight0, n, grids[0]);
+
+ GlobalCoords lowerLeft1({0, 1});
+ GlobalCoords upperRight1({2, 2});
+ buildGrid(lowerLeft1, upperRight1, n, grids[1]);
+
+ // writing grids
+ for (size_t i=0; i<grids.size(); i++) {
+ const auto& levelView = grids[i]->levelGridView(0);
+ writeToVTK(levelView, path, "body_" + std::to_string(i) + "_level0");
+ }
+
+ // compute coupling boundaries
+ LevelView gridView0 = grids[0]->levelGridView(0);
+ LevelView gridView1 = grids[1]->levelGridView(0);
+ LevelBoundaryPatch upper(gridView0);
+ LevelBoundaryPatch lower(gridView1);
+
+ lower.insertFacesByProperty([&](typename LevelView::Intersection const &in) {
+ return xyBetween(lowerLeft1, {upperRight1[0], lowerLeft1[1]}, in.geometry().center());
+ });
+
+ upper.insertFacesByProperty([&](typename LevelView::Intersection const &in) {
+ return xyBetween({lowerLeft0[0], upperRight0[1]}, upperRight0, in.geometry().center());
+ });
+
+ // set contact coupling
+ Dune::Contact::NormalProjection<LeafBoundaryPatch> contactProjection;
+
+ CouplingPair coupling;
+ coupling.set(0, 1, upper, lower, 0.1, CouplingPair::CouplingType::STICK_SLIP, contactProjection, nullptr);
+
+ double coveredArea_ = 0.8;
+ CouplingType contactCoupling;
+ contactCoupling.setGrids(*grids[0], *grids[1]);
+ contactCoupling.setupContactPatch(*coupling.patch0(),*coupling.patch1());
+ contactCoupling.gridGlueBackend_ = coupling.backend();
+ contactCoupling.setCoveredArea(coveredArea_);
+ contactCoupling.setup();
+
+ using Basis = P1NodalBasis<LevelView, c_type>;
+ Basis basis0(grids[0]->levelGridView(0));
+ Basis basis1(grids[1]->levelGridView(0));
+
+ printBasisDofLocation(basis0);
+ printBasisDofLocation(basis1);
+
+ ScalarBlockVector nWeights0, nWeights1;
+ NodalWeights<Basis, Basis> nodalWeights(basis0, basis1);
+ nodalWeights.assemble(*contactCoupling.getGlue(), nWeights0, nWeights1, true);
+
+ ScalarBlockVector weights0;
+ {
+ Assembler<Basis, Basis> assembler(basis0, basis0);
+ NeumannBoundaryAssembler<GridType, typename ScalarBlockVector::block_type> boundaryAssembler(
+ std::make_shared<ConstantFunction<
+ GlobalCoords, typename ScalarBlockVector::block_type>>(1));
+ assembler.assembleBoundaryFunctional(boundaryAssembler, weights0, upper);
+ }
+
+ ScalarBlockVector weights1;
+ {
+ Assembler<Basis, Basis> assembler(basis1, basis1);
+ NeumannBoundaryAssembler<GridType, typename ScalarBlockVector::block_type> boundaryAssembler(
+ std::make_shared<ConstantFunction<
+ GlobalCoords, typename ScalarBlockVector::block_type>>(1));
+ assembler.assembleBoundaryFunctional(boundaryAssembler, weights1, lower);
+ }
+
+ print(weights0, "assembled weights0: ");
+ print(nWeights0, "nWeights0: ");
+
+ print(weights1, "assembled weights1: ");
+ print(nWeights1, "nWeights1: ");
+
+ bool sizePassed = true;
+ if (weights0.size() != nWeights0.size() | weights1.size() != nWeights1.size()) {
+ sizePassed = false;
+ }
+
+ bool entriesPassed0 = true;
+ for (size_t i=0; i<weights0.size(); i++) {
+ entriesPassed0 = entriesPassed0 & isClose(weights0[i], nWeights0[i]);
+ }
+
+ bool entriesPassed1 = true;
+ for (size_t i=0; i<weights1.size(); i++) {
+ entriesPassed1 = entriesPassed1 & isClose(weights1[i], nWeights1[i]);
+ }
+
+ bool passed = sizePassed & entriesPassed0 & entriesPassed1;
+
+ std::cout << "Overall, the test " << (passed ? "was successful!" : "failed!") << std::endl;
+
+ std::cout.rdbuf(coutbuf); //reset to standard output again
+ return passed ? 0 : 1;
+
+} catch (Dune::Exception &e) {
+ Dune::derr << "Dune reported error: " << e << std::endl;
+} catch (std::exception &e) {
+ std::cerr << "Standard exception: " << e.what() << std::endl;
+} // end try
+} // end main
diff --git a/src/utils/debugutils.hh b/src/utils/debugutils.hh
index 150aad1c..c81de410 100644
--- a/src/utils/debugutils.hh
+++ b/src/utils/debugutils.hh
@@ -12,7 +12,7 @@
#include <dune/grid/io/file/vtk/vtkwriter.hh>
-using namespace std;
+//using namespace std;
template <int s>
@@ -33,7 +33,7 @@ using namespace std;
std::cout << message << std::endl;
for (size_t i=0; i<x.size(); i++) {
- std::cout << x[i] << ";" << std::endl;;
+ std::cout << x[i] << ";" << std::endl;
}
}
@@ -184,9 +184,45 @@ using namespace std;
std::cout.rdbuf( lBufferOld );
}
+ template <class Intersection, class Basis>
+ void printIntersection(const Intersection& it, const Basis& basis, std::string message = "") {
+ if (message != "") {
+ std::cout << message << std::endl;
+ }
+
+ const auto& inside = it.inside();
+ const auto& localCoefficients = basis.getLocalFiniteElement(inside).localCoefficients();
+
+ std::cout << "dofs: ";
+ for (size_t i=0; i<localCoefficients.size(); i++) {
+ unsigned int entity = localCoefficients.localKey(i).subEntity();
+ unsigned int codim = localCoefficients.localKey(i).codim();
+
+ if (it.containsInsideSubentity(entity, codim)) {
+ int dofIndex = basis.index(it.inside(), i);
+ std::cout << dofIndex << " ";
+ break;
+ }
+ }
+ std::cout << std::endl;
+ }
+
+ template <class BoundaryPatch, class Basis>
+ void printBoundary(const BoundaryPatch& patch, const Basis& basis, std::string message = "") {
+ if (message != "") {
+ std::cout << message << std::endl;
+ }
+
+ for (auto bIt = patch.begin(); bIt != patch.end(); ++bIt) {
+ printIntersection(*bIt, basis, "");
+ }
+
+ std::cout << std::endl;
+ }
+
template <class BasisType, typename ctype=double>
void printBasisDofLocation(const BasisType& basis) {
- typedef typename BasisType::GridView GridView;
+ /* typedef typename BasisType::GridView GridView;
const int dim = GridView::dimension;
@@ -197,12 +233,9 @@ using namespace std;
const GridView& gridView = basis.getGridView();
const int gridN = std::pow(gridView.indexSet().size(dim), 1.0/dim)-1;
- typedef typename GridView::template Codim<0>::Iterator ElementIterator;
- ElementIterator it = gridView.template begin<0>();
- ElementIterator end = gridView.template end<0>();
- for (; it != end; ++it) {
- const typename BasisType::LocalFiniteElement& tFE = basis.getLocalFiniteElement(*it);
- const auto& geometry = (*it).geometry();
+ for (auto& it : elements(gridView)) {
+ const typename BasisType::LocalFiniteElement& tFE = basis.getLocalFiniteElement(it);
+ const auto& geometry = it.geometry();
for(int i=0; i<geometry.corners(); ++i) {
const auto& vertex = geometry.corner(i);
@@ -219,7 +252,7 @@ using namespace std;
for(int k=0; k<localBasisSize; ++k) {
if (localBasisRep[k]==1) {
- int dofIndex = basis.index((*it), k);
+ int dofIndex = basis.index(it, k);
if (indexTransformation.count(dofIndex)==0) {
indexTransformation[dofIndex] = vertexIndex;
@@ -245,5 +278,50 @@ using namespace std;
}
std::cout << std::endl;
}
+ std::cout << std::endl;*/
+
+ const int dim = BasisType::GridView::dimension;
+
+ std::map<int, Dune::FieldVector<ctype, dim>> indexCoords;
+
+
+ const auto& gridView = basis.getGridView();
+
+ for (auto& it : elements(gridView)) {
+ const typename BasisType::LocalFiniteElement& tFE = basis.getLocalFiniteElement(it);
+ const auto& geometry = it.geometry();
+
+ for(int i=0; i<geometry.corners(); ++i) {
+ const auto& vertex = geometry.corner(i);
+ const auto& local = geometry.local(vertex);
+
+ const int localBasisSize = tFE.localBasis().size();
+ std::vector<Dune::FieldVector<double, 1>, std::allocator<Dune::FieldVector<double, 1> > > localBasisRep(localBasisSize);
+ tFE.localBasis().evaluateFunction(local, localBasisRep);
+
+ for(int k=0; k<localBasisSize; ++k) {
+ if (1.0 - localBasisRep[k]< 1e-14) {
+ int dofIndex = basis.index(it, k);
+ indexCoords[dofIndex] = vertex;
+ break;
+ }
+ }
+ }
+ }
+
+ std::cout << "Coords of basis dofs: " << indexCoords.size() << std::endl;
+ auto&& mapIt = indexCoords.begin();
+ const auto& mapEndIt = indexCoords.end();
+ for (; mapIt!=mapEndIt; ++mapIt) {
+ std::cout << mapIt->first << " Coords: ";
+
+ const auto& coords = mapIt->second;
+ for (size_t i=0; i<coords.size(); i++) {
+ std::cout << coords[i] << " ";
+ }
+ std::cout << std::endl;
+ }
+
+ std::cout << std::endl;
}
#endif
--
GitLab