From 0f474e21931493a3dc8e3b6f99b95542fb2b3a22 Mon Sep 17 00:00:00 2001
From: Klaus Boehnlein <klaus.boehnlein@tu-dresden.de>
Date: Mon, 12 Jul 2021 17:31:11 +0200
Subject: [PATCH] Update for Namespace Experimental/PeriodicBasis
---
dune/microstructure/prestrainimp.hh | 229 ++++++++++------------------
src/dune-microstructure.cc | 4 +-
2 files changed, 83 insertions(+), 150 deletions(-)
diff --git a/dune/microstructure/prestrainimp.hh b/dune/microstructure/prestrainimp.hh
index e5268c68..a3b68f0a 100644
--- a/dune/microstructure/prestrainimp.hh
+++ b/dune/microstructure/prestrainimp.hh
@@ -1,29 +1,42 @@
-#ifndef DUNE_REDUCEDROD_PRESTRAINIMP_HH
-#define DUNE_REDUCEDROD_PRESTRAINIMP_HH
+#ifndef DUNE_MICROSTRUCTURE_PRESTRAINIMP_HH
+#define DUNE_MICROSTRUCTURE_PRESTRAINIMP_HH
-#include <dune/microstructure/microstructuredrodgrid.hh>
+// #include <dune/microstructure/microstructuredrodgrid.hh>
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/yaspgrid.hh>
//template<class MicrostructuredRodGrid> reicht include
+
+using namespace Dune;
+
+
class PrestrainImp {
public:
+ static const int dim = 3;
static const int nCompo = 3;
- using GridType_Ce = typename MicrostructuredRodGrid::CellGridType;
- using Domain = typename GridType_Ce::LeafGridView::template Codim<0>::Geometry::GlobalCoordinate;
- using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
- using Func2Tensor = std::function< MatrixRT(const Domain&) >;
+// using GridType_Ce = typename MicrostructuredRodGrid::CellGridType;
+ using CellGridType = YaspGrid< dim, EquidistantOffsetCoordinates< double, dim>>;
+ using GridView = CellGridType::LeafGridView;
+ using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
+
+ using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
+ using Func2Tensor = std::function< MatrixRT(const Domain&) >;
+// using Domain = typename GridType_Ce::LeafGridView::template Codim<0>::Geometry::GlobalCoordinate;
+// using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
+// using Func2Tensor = std::function< MatrixRT(const Domain&) >;
protected:
- double p1, p2, a;
+ double p1, p2, theta;
double width; //cell geometry
public:
- PrestrainImp(double p1, double p2, double a, double width) : p1(p1), p2(p2), a(a), width(width){}
+ PrestrainImp(double p1, double p2, double theta, double width) : p1(p1), p2(p2), theta(theta), width(width){}
Func2Tensor getPrestrain(unsigned int imp)
{
@@ -33,149 +46,69 @@ public:
using std::sin;
using std::cos;
- if (imp==1) // Prestrain inducing twist (bisher nur extension)
- {
-
- Func2Tensor B1_ = [this] (const Domain& z)
- {
- MatrixRT ret = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}};
- if (z[0] < 0.25 && z[1] < 0 ||
- 0.25 < z[0] && z[0] < 0.5 && z[2] < 0 ||
- 0.5 < z[0] && z[0] < 0.75 && z[1] >= 0 ||
- 0.75 < z[0] && z[2] >= 0)
- ret *= this->p1;
- else
- ret *= this->p2;
- return ret;
- };
- return B1_;
- }
- else if (imp==2) // Splay bend. B = 1/6*a*Id + 1/2*a*(n otimes n). n= (cos, 0, sin)
- {
-
- Func2Tensor B2_ = [this] (const Domain& z)
- {
- auto pi = std::acos(-1.0);
- double beta = pi/4.0 + pi/4.0*z[1]; //z[1]=x3
- MatrixRT Id(0);
- for (int i=0;i<nCompo;i++)
- Id[i][i]=1.0;
- MatrixRT pressure = Id;
- pressure *= 1.0/6.0;
- MatrixRT n_ot_n = {
- {cos(beta)*cos(beta), 0.0, cos(beta)*sin(beta)},
- {0.0, 0.0, 0.0 },
- {cos(beta)*sin(beta), 0.0, sin(beta)*sin(beta)}
- };
- n_ot_n*=0.5;
- pressure += n_ot_n;
- pressure *= this->a;
- return pressure;
- };
- return B2_;
- }
- else if (imp==3) // Twist
- {
-
- Func2Tensor B3_ = [this] (const Domain& z)
- {
- auto pi = std::acos(-1.0);
- double beta = pi/4.0 + pi/4.0*z[1];
- MatrixRT pressure = {{1.0/6.0, 0.0, 0.0}, {0.0, 1.0/6.0, 0.0}, {0.0, 0.0, 1.0/6.0}};
- MatrixRT n_ot_n = {
- {0.5 * cos(beta)*cos(beta), 0.5 * cos(beta)*sin(beta), 0.0},
- {0.5 * cos(beta)*sin(beta), 0.5 * sin(beta)*sin(beta), 0.0},
- {0.0, 0.0, 0.0 }
- };
- pressure += n_ot_n;
- pressure *= this->a;
- return pressure;
- };
- return B3_;
- }
- else if (imp==4) // Splay bend different directions
- {
-
- Func2Tensor B4_ = [this] (const Domain& z)
- {
- auto pi = std::acos(-1.0);
- MatrixRT pressure = {{1.0/6.0, 0.0, 0.0}, {0.0, 1.0/6.0, 0.0}, {0.0, 0.0, 1.0/6.0}};
- if ( z[2] < 0.5 ){
- double beta = pi/4.0 + pi/4.0*z[1];
- MatrixRT n_ot_n = {
- {0.5 * cos(beta)*cos(beta), 0.0, 0.5 * cos(beta)*sin(beta)}, {0.0, 0.0, 0.0}, {0.5 * cos(beta)*sin(beta), 0.0, 0.5* sin(beta)*sin(beta)}
- };
- pressure += n_ot_n;
- }
- else{
- double beta = pi/4.0 + pi/4.0*z[0];
- MatrixRT n_ot_n = {
- {0.5 * cos(beta)*cos(beta), 0.5 * cos(beta)*sin(beta), 0.0}, {0.5 * cos(beta)*sin(beta), 0.5* sin(beta)*sin(beta), 0.0}, {0.0, 0.0, 0.0}
- };
- pressure += n_ot_n;
- }
- pressure *= this->a;
- return pressure;
- };
- return B4_;
- }
- else if (imp==5) // Mix twist and splay bend
- {
-
- Func2Tensor B5_ = [this] (const Domain& z)
- {
- auto pi = std::acos(-1.0);
- double beta = pi/4.0 + pi/4.0*z[1];
- MatrixRT pressure = {{1.0/6.0, 0.0, 0.0}, {0.0, 1.0/6.0, 0.0}, {0.0, 0.0, 1.0/6.0}};
- if ( z[2] < 0.5 ){
- MatrixRT n_ot_n = {
- {0.5 * cos(beta)*cos(beta), 0.0, 0.5 * cos(beta)*sin(beta)}, {0.0, 0.0, 0.0}, {0.5 * cos(beta)*sin(beta), 0.0, 0.5* sin(beta)*sin(beta)}
- };
- pressure += n_ot_n;
- }
- else{
- MatrixRT n_ot_n = {
- {0.5 * cos(beta)*cos(beta), 0.5 * cos(beta)*sin(beta), 0.0}, {0.5 * cos(beta)*sin(beta), 0.5* sin(beta)*sin(beta), 0.0}, {0.0, 0.0, 0.0}
- };
- pressure += n_ot_n;
- }
- pressure *= this->a;
- return pressure;
- };
- return B5_;
- }
- else if (imp==6) // Shear prestrain
+ if (imp==1)
{
- Func2Tensor B6_ = [this](const Domain& z) {
- MatrixRT ret = {
- {0.0, 0.5, 0.0},
- {0.5, 0.0, 0.0},
- {0.0, 0.0, 0.0}};
- ret *= this->a;
- return ret;
- };
- return B6_;
- }
- else // Prestrain inducing twist (ist ein Versuch)
- {
-
- Func2Tensor B7_ = [this] (const Domain& z)
- {
- MatrixRT ret = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}};
- if (z[0] < 0.5 && sqrt(pow(z[1],2) + pow(z[2],2) ) < this->width/4.0 ||
- 0.5 < z[0] && sqrt(pow(z[1],2) + pow(z[2],2) ) > this->width/4.0)
- ret *= this->p1;
- else
- ret *= this->p2;
- return ret;
- };
- return B7_;
- }
-
+ Func2Tensor B1_ = [this] (const Domain& x) { // ISOTROPIC PRESSURE
+ if (abs(x[0]) > (theta/2) && x[2] > 0)
+ return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+ if (abs(x[0]) < (theta/2) && x[2] < 0)
+ return MatrixRT{{p2, 0.0 , 0.0}, {0.0, p2, 0.0}, {0.0, 0.0, p2}};
+ else
+ return MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+
+ };
+ std::cout <<" Prestrain Type: 1 "<< std::endl;
+ return B1_;
+ }
+ else if (imp==2)
+ {
+ Func2Tensor B2_ = [this] (const Domain& x) { // ISOTROPIC PRESSURE
+
+ if (abs(x[0]) < (theta/2) && x[2] < 0 && x[2] > -(1.0/2.0) )
+ return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+ else
+ return MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+
+
+ };
+ std::cout <<" Prestrain Type: 2 "<< std::endl;
+ return B2_;
+
+
+ }
+
+
+
+
+ // TODO ANISOTROPIC PRESSURE
+
+// else if (imp==2)
+// {
+//
+// Func2Tensor B2_ = [this] (const Domain& z)
+// {
+// auto pi = std::acos(-1.0);
+// double beta = pi/4.0 + pi/4.0*z[1]; //z[1]=x3
+// MatrixRT Id(0);
+// for (int i=0;i<nCompo;i++)
+// Id[i][i]=1.0;
+// MatrixRT pressure = Id;
+// pressure *= 1.0/6.0;
+// MatrixRT n_ot_n = {
+// {cos(beta)*cos(beta), 0.0, cos(beta)*sin(beta)},
+// {0.0, 0.0, 0.0 },
+// {cos(beta)*sin(beta), 0.0, sin(beta)*sin(beta)}
+// };
+// n_ot_n*=0.5;
+// pressure += n_ot_n;
+// pressure *= this->a;
+// return pressure;
+// };
+// return B2_;
+// }
}
-
};
diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index f7d609b7..c49ad866 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -1394,7 +1394,7 @@ int main(int argc, char *argv[])
/////////////////////////////////////////////////////////
using namespace Functions::BasisFactory;
- PeriodicIndexSet periodicIndices;
+ Functions::Experimental::BasisFactory::PeriodicIndexSet periodicIndices;
// Don't do the following in real life: It has quadratic run-time in the number of vertices.
for (const auto& v1 : vertices(gridView_CE))
@@ -1406,7 +1406,7 @@ int main(int argc, char *argv[])
gridView_CE,
power<dim>(
// lagrange<1>(),
- periodic(lagrange<1>(), periodicIndices),
+ Functions::Experimental::BasisFactory::periodic(lagrange<1>(), periodicIndices),
flatLexicographic()));//
// flatInterleaved()));
--
GitLab