Loading [MathJax]/extensions/tex2jax.js

dune-composites (unstable)

linearStaticDriver.hh
1// -*- tab-width: 4; indent-tabs-mode: nil -*-
2// vi: set et ts=4 sw=4 sts=4:
3#include <dune/composites/Setup/parallelPartition.hh>
4#include <dune/composites/Setup/GridTransformation.hh>
5#include <dune/composites/Setup/dirichletBC.hh>
6#include <dune/composites/Driver/FEM/Serendipity/serendipityfem.hh>
7
8#include "../PostProcessing/computeStresses.hh"
9#include "../PostProcessing/plot_properties.hh"
10#include "localOperator/linearelasticity.hh"
11
12namespace Dune{
13 namespace Composites{
14
16 template <typename MODEL>
17 class linearStaticDriver{
18
19 public:
20 typedef Dune::PDELab::ISTL::BCRSMatrixBackend<> MBE;
21 typedef double RF;
22
23 int elem_order;
24
25 Dune::MPIHelper& helper;
26
27 linearStaticDriver(Dune::MPIHelper & helper_) : helper(helper_){
28 elem_order = 2; // Default Value which Can be overwritten by setElementOrder
29 };
30
31 void inline setElementOrder(int val = 2){
32 elem_order = val;
33 }
34
35 void inline apply(MODEL& myModel){
36 if (myModel.setUp_Required == true){
37 myModel.LayerCake();
38 myModel.setUpMaterials();
39 myModel.computeElasticTensors();
40 }
41
42 Dune::Timer watch;
43 std::vector<double> times(3);
44
45 // === A: Setup YaspGrid
46 watch.reset();
47 typedef Dune::YaspGrid<3,Dune::TensorProductCoordinates<RF,3>> YGRID;
48 YaspPartition<3> yp;
49 YGRID yaspgrid(myModel.coords,myModel.periodic,myModel.overlap,helper.getCommunicator(),&yp);
50 int refinements = myModel.refineBaseGrid;
51 if(refinements > 0) yaspgrid.globalRefine(refinements);
52 typedef YGRID::LeafGridView YGV;
53 const YGV ygv = yaspgrid.leafGridView();
54
55 if(helper.rank() == 0){
56 std::cout << "Number of elements per processor: " << ygv.size(0) << std::endl;
57 std::cout << "Number of nodes per processor: " << ygv.size(3) << std::endl;
58 }
59
60 myModel.setPG(ygv); // Loops over elements and assigns a physical group
61
62 // ==================================================================
63 // Transform Grid
64 // ==================================================================
65 typedef Dune::Composites::GridTransformation<3,MODEL> GRID_TRAFO;
66 GRID_TRAFO gTrafo(myModel, helper.rank());
67
68 typedef typename Dune::GeometryGrid<YGRID,GRID_TRAFO> GRID;
69 GRID grid(yaspgrid,gTrafo);
70 if(helper.rank() == 0)
71 std::cout << "Grid transformation complete" << std::endl;
72
73 //Define Grid view
74 typedef typename GRID::LeafGridView GV;
75 const GV gv = grid.leafGridView();
76
77 if(helper.rank() == 0)
78 std::cout << "Grid view set up" << std::endl;
79
80 times[1] = watch.elapsed();
81 watch.reset();
82
83 typedef Scalar_BC<GV,MODEL,RF> BC;
84 typedef Dune::PDELab::CompositeConstraintsParameters<BC,BC,BC> Constraints;
85
86 // ==================================================================
87 // Set up problem
88 // ==================================================================
89 myModel.template computeTensorsOnGrid<GV,YGV>(gv,ygv);
90
91 // Setup initial boundary conditions for each degree of freedom
92 typedef Scalar_BC<GV,MODEL,RF> InitialDisp;
93 InitialDisp u1(gv, myModel,0), u2(gv, myModel,1), u3(gv, myModel,2);
94
95 // Wrap scalar boundary conditions in to vector
96 typedef Dune::PDELab::CompositeGridFunction<InitialDisp,InitialDisp,InitialDisp> InitialSolution;
97 InitialSolution initial_solution(u1,u2,u3);
98
99 // Construct grid function spaces for each degree of freedom
100 typedef Dune::PDELab::OverlappingConformingDirichletConstraints CON;
101 CON con;
102 typedef Dune::PDELab::ConformingDirichletConstraints CON_EXT; //needed only for Geneo solver
103
104 times[2] = watch.elapsed();
105 watch.reset();
106
107 if(elem_order == 2){ // == Element Order
108 const int element_order = 2;
109 const int dofel = 3 * 20;
110 const int non_zeros = 81;
111
112 if(helper.rank() == 0)
113 std::cout << "Piecewise quadratic serendipity elements" << std::endl;
114
115 typedef Dune::PDELab::SerendipityLocalFiniteElementMap<GV,typename GV::Grid::ctype,RF,element_order> FEM;
116 FEM fem(gv);
117
118 // Create AllEntitySet to be passed to ALL grid function spaces
119 typedef Dune::PDELab::AllEntitySet<GV> ES_ALL;
120 ES_ALL es_all(gv);
121
122 using SCALAR_VBE = Dune::PDELab::ISTL::VectorBackend<>;
123 typedef Dune::PDELab::GridFunctionSpace<ES_ALL, FEM, CON, SCALAR_VBE> SCALAR_GFS;
124 typedef Dune::PDELab::GridFunctionSpace<ES_ALL, FEM, CON_EXT, SCALAR_VBE> SCALAR_GFS_EXT;
125 SCALAR_GFS U1(es_all,fem,con); U1.name("U1"); SCALAR_GFS_EXT U1_EXT(es_all,fem);
126 SCALAR_GFS U2(es_all,fem,con); U2.name("U2"); SCALAR_GFS_EXT U2_EXT(es_all,fem);
127 SCALAR_GFS U3(es_all,fem,con); U3.name("U3"); SCALAR_GFS_EXT U3_EXT(es_all,fem);
128
129 // Note that Vectors are blocked by Dune::PDELab::EntityBlockedOrderingTag
130 const int block_size = 3;
131 using VBE = Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed>;
132 using Ordering = Dune::PDELab::EntityBlockedOrderingTag;
133 typedef Dune::PDELab::PowerGridFunctionSpace <SCALAR_GFS,block_size,VBE,Ordering> GFS;
134 const GFS gfs(U1,U2,U3);
135 typedef Dune::PDELab::PowerGridFunctionSpace <SCALAR_GFS_EXT,block_size,VBE,Ordering> GFS_EXT;
136 const GFS_EXT gfs_ext(U1_EXT,U2_EXT,U3_EXT);
137
138 typedef typename GFS::template ConstraintsContainer<RF>::Type C;
139
140 // Make constraints map and initialize it from a function
141 C cg;
142 cg.clear();
143
144 BC U1_cc(gv,myModel,0), U2_cc(gv,myModel,1), U3_cc(gv,myModel,2);
145 Constraints constraints(U1_cc,U2_cc,U3_cc);
146
147 Dune::PDELab::constraints(constraints,gfs,cg);
148
149 MBE mbe(non_zeros); // Maximal number of nonzeroes per row
150
151 // === Construct Linear Operator on FEM Space
152 typedef Dune::PDELab::linearelasticity<GV, MODEL, dofel> LOP;
153 LOP lop(gv, myModel);
154
155 typedef Dune::PDELab::GridOperator<GFS,GFS,LOP,MBE,RF,RF,RF,C,C> GO;
156 GO go(gfs,cg,gfs,cg,lop,mbe);
157 typedef Dune::PDELab::GridOperator<GFS_EXT,GFS_EXT,LOP,MBE,RF,RF,RF,C,C> GO_EXT;
158
159 // === Make coefficent vector and initialize it from a function
160 typedef Dune::PDELab::Backend::Vector<GFS,double> V;
161 V u(gfs,0.0);
162 Dune::PDELab::interpolate(initial_solution,gfs,u);
163
164 // === Set non constrained dofs to zero
165 Dune::PDELab::set_shifted_dofs(cg,0.0,u);
166 myModel.template solve<GO,GO_EXT,V,GFS,GFS_EXT,C,Constraints,MBE,LOP>(go,u,gfs,gfs_ext,cg,constraints,mbe,lop);
167
168 Dune::SubsamplingVTKWriter<GV> vtkwriter(gv,Dune::refinementLevels(0));
169 Dune::PDELab::addSolutionToVTKWriter(vtkwriter,gfs,u);
170 vtkwriter.write(myModel.vtk_displacement_output,Dune::VTK::appendedraw);
171
172 calculateStresses<MODEL,V,GV,GFS,MBE>(myModel,u,gv,gfs,mbe);
173
174 plotProperties<MODEL,V,GV,GFS,MBE>(myModel,gv,gfs,mbe);
175
176 myModel.template postprocess<GO,V,GFS,C,MBE,GV>(go,u,gfs,cg,gv,mbe);
177 }
178 } // end apply
179
180 };
181 }
182}
Dune::Composites::GridTransformation
Grid Transformation.
Definition: GridTransformation.hh:19
Dune::Composites::Scalar_BC
Define Scalar Dirichlet Boundary Conditions.
Definition: dirichletBC.hh:14
Dune::Composites::YaspPartition
Partition yaspgrid for parallelism.
Definition: parallelPartition.hh:12
Dune::Composites::linearStaticDriver
linear static driver
Definition: linearStaticDriver.hh:17
Dune::PDELab::SerendipityLocalFiniteElementMap
Definition: serendipityfem.hh:20
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden & Uni Heidelberg  |  generated with Hugo v0.111.3 (Apr 3, 22:46, 2025)