dune-composites (2.5.1)

heatStaticDriver.hh
1#include <dune/composites/Setup/parallelPartition.hh>
2#include <dune/composites/Setup/GridTransformation.hh>
3#include <dune/composites/Setup/dirichletBC.hh>
4#include <dune/composites/Driver/FEM/Serendipity/serendipityfem.hh>
5
6#include "../PostProcessing/computeStresses.hh"
7#include "../PostProcessing/plot_properties.hh"
8#include "localOperator/linearelasticity.hh"
9
10namespace Dune{
11 namespace Composites{
12
14
18 template <typename MODEL>
19 class heatStaticDriver{
20
21 public:
22 typedef Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::none,1> Scalar_VectorBackend;
23 typedef Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::fixed,3> VectorBackend;
24 typedef Dune::PDELab::istl::BCRSMatrixBackend<> MBE;
25 typedef double RF;
26
27 int elem_order;
28
29 Dune::MPIHelper& helper;
30
31 heatStaticDriver(Dune::MPIHelper & helper_) : helper(helper_){
32 elem_order = 2; // Default Value which Can be overwritten by setElementOrder
33 };
34
35 void inline setElementOrder(int val = 2){
36 elem_order = val;
37 }
38
39 void inline apply(MODEL& myModel){
40 if (myModel.setUp_Required == true){
41 myModel.LayerCake();
42 myModel.setUpMaterials();
43 myModel.computeElasticTensors();
44 }
45
46 Dune::Timer watch;
47 std::vector<double> times(3);
48
49 // === A: Setup YaspGrid
50 watch.reset();
51 typedef Dune::YaspGrid<3,Dune::TensorProductCoordinates<double,3>> YGRID;
52 YaspPartition<3> yp(helper.size(),myModel.GridPartition());
53 YGRID yaspgrid(myModel.coords,myModel.periodic,myModel.overlap,helper.getCommunicator(),(Dune::YLoadBalance<3>*)&yp);
54 int refinements = myModel.refineBaseGrid;
55 if(refinements > 0) yaspgrid.globalRefine(refinements);
56 typedef YGRID::LeafGridView YGV;
57 const YGV& ygv = yaspgrid.leafGridView();
58 int size = yaspgrid.globalSize(0)*yaspgrid.globalSize(1)*yaspgrid.globalSize(2);
59
60 if(helper.rank() == 0){
61 std::cout << "Number of elements per processor: " << ygv.size(0) << std::endl;
62 std::cout << "Number of nodes per processor: " << ygv.size(3) << std::endl;
63 }
64
65 myModel.setPG(ygv); // Loops over elements and assigns a physical group
66
67 // ==================================================================
68 // Transform Grid
69 // ==================================================================
70 typedef Dune::Composites::GridTransformation<3,MODEL> GRID_TRAFO;
71 GRID_TRAFO gTrafo(myModel,helper.rank());
72
73 typedef typename Dune::GeometryGrid<YGRID,GRID_TRAFO> GRID;
74 GRID grid(yaspgrid,gTrafo);
75 if(helper.rank() == 0)
76 std::cout << "Grid transformation complete" << std::endl;
77
78 //Define Grid view
79 typedef typename GRID::LeafGridView GV;
80 const GV& gv = grid.leafGridView();
81
82 if(helper.rank() == 0)
83 std::cout << "Grid view set up" << std::endl;
84
85 times[1] = watch.elapsed();
86 watch.reset();
87
88 typedef Scalar_BC<GV,MODEL,RF> BC;
89 typedef Dune::PDELab::CompositeConstraintsParameters<BC,BC,BC> Constraints;
90
91 // ==================================================================
92 // Set up problem
93 // ==================================================================
94 myModel.template computeTensorsOnGrid<GV,YGV>(gv,ygv);
95
96 // Setup initial boundary conditions for each degree of freedom
97 typedef Scalar_BC<GV,MODEL,double> InitialDisp;
98 InitialDisp u1(gv, myModel,0), u2(gv, myModel,1), u3(gv, myModel,2);
99
100 // Wrap scalar boundary conditions in to vector
101 typedef Dune::PDELab::CompositeGridFunction<InitialDisp,InitialDisp,InitialDisp> InitialSolution;
102 InitialSolution initial_solution(u1,u2,u3);
103
104 // Construct grid function spaces for each degree of freedom
105 typedef Dune::PDELab::OverlappingConformingDirichletConstraints CON;
106 CON con;
107
108 times[2] = watch.elapsed();
109 watch.reset();
110
111 if(elem_order == 2){ // == Element Order
112 const int element_order = 2;
113 const int dofel = 3 * 20;
114 const int non_zeros = 81;
115
116 if(helper.rank() == 0)
117 std::cout << "Piecewise quadratic serendipity elements" << std::endl;
118
119 typedef Dune::PDELab::SerendipityLocalFiniteElementMap<GV,typename GV::Grid::ctype,double,element_order> FEM;
120 FEM fem(gv);
121
122 typedef Dune::PDELab::GridFunctionSpace<GV, FEM, CON, Scalar_VectorBackend> SCALAR_GFS;
123 SCALAR_GFS dispU1(gv,fem,con); dispU1.name("U1");
124 SCALAR_GFS dispU2(gv,fem,con); dispU2.name("U2");
125 SCALAR_GFS dispU3(gv,fem,con); dispU3.name("U3");
126
127 // Note that Vectors are blocked by Dune::PDELab::EntityBlockedOrderingTag
128 typedef Dune::PDELab::CompositeGridFunctionSpace <VectorBackend,Dune::PDELab::EntityBlockedOrderingTag, SCALAR_GFS, SCALAR_GFS, SCALAR_GFS> GFS;
129 const GFS gfs(dispU1, dispU2, dispU3);
130
131 typedef typename GFS::template ConstraintsContainer<RF>::Type C;
132
133 // Make constraints map and initialize it from a function
134 C cg;
135 cg.clear();
136
137 BC U1_cc(gv,myModel,0), U2_cc(gv,myModel,1), U3_cc(gv,myModel,2);
138 Constraints constraints(U1_cc,U2_cc,U3_cc);
139 Dune::PDELab::constraints(constraints,gfs,cg);
140
141 MBE mbe(non_zeros); // Maximal number of nonzeroes per row
142
143 //TODO precompute geneo eigenvectors and reuse
144 //Initial idea, compute pure thermal, pure mechanical to calculate a good initial guess for the failure
145 //load. Turns out at least in this simple case to require more solves than just guessing and iterating
146 // === Construct Linear Operator on FEM Space
147 // 1. Pure thermal
148 /*double therm_fail = 0;
149 {
150 myModel.setThermal(true);
151 myModel.setPressure(0.0);
152
153 typedef Dune::PDELab::linearelasticity<GV, MODEL, dofel> LOP;
154 LOP lop(gv, myModel);
155
156 typedef Dune::PDELab::GridOperator<GFS,GFS,LOP,MBE,RF,RF,RF,C,C> GO;
157 GO go(gfs,cg,gfs,cg,lop,mbe);
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
167 // === Solve the linear system
168 myModel.template solve<GO,V,GFS,C,Constraints,MBE,LOP>(go,u,gfs,cg,constraints,mbe,lop);
169
170 if(gv.comm().rank() == 0) std::cout << "=== Calculate Stresses" << std::endl;
171 calculateStresses<MODEL,V,GV,GFS,MBE>(myModel,u,gv,gfs,mbe);
172
173 // Pure thermal failure
174 using Dune::PDELab::Backend::native;
175 double max = 0.0;
176 typedef typename GV::Traits::template Codim<0>::Iterator ElementIterator;
177 for (ElementIterator it = gv.template begin<0>(); it!=gv.template end<0>(); ++it)
178 { // loop through each element
179 int id = gv.indexSet().index(*it);
180 auto stress = myModel.getStress(id);
181 double h = myModel.FailureCriteria(stress);
182 if(h > max) max = h;
183 }
184 MPI_Allreduce(&max, &therm_fail, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
185 if(gv.comm().rank() == 0) std::cout << "Pure thermal failure " << therm_fail << std::endl;
186 }
187 //Now get the maximum failure criterion here
188 // 2. Run without thermal load but with
189 double mech_fail = 0;
190 {
191 myModel.setThermal(false);
192 myModel.setPressure(1.0);
193
194 typedef Dune::PDELab::linearelasticity<GV, MODEL, dofel> LOP;
195 LOP lop(gv, myModel);
196
197 typedef Dune::PDELab::GridOperator<GFS,GFS,LOP,MBE,RF,RF,RF,C,C> GO;
198 GO go(gfs,cg,gfs,cg,lop,mbe);
199
200 // === Make coefficent vector and initialize it from a function
201 typedef Dune::PDELab::Backend::Vector<GFS,double> V;
202 V u(gfs,0.0);
203 Dune::PDELab::interpolate(initial_solution,gfs,u);
204
205 // === Set non constrained dofs to zero
206 Dune::PDELab::set_shifted_dofs(cg,0.0,u);
207
208 // === Solve the linear system
209 myModel.template solve<GO,V,GFS,C,Constraints,MBE,LOP>(go,u,gfs,cg,constraints,mbe,lop);
210
211 if(gv.comm().rank() == 0) std::cout << "=== Calculate Stresses" << std::endl;
212 calculateStresses<MODEL,V,GV,GFS,MBE>(myModel,u,gv,gfs,mbe);
213
214 // Pure thermal failure
215 using Dune::PDELab::Backend::native;
216 double max = 0.0;
217 typedef typename GV::Traits::template Codim<0>::Iterator ElementIterator;
218 for (ElementIterator it = gv.template begin<0>(); it!=gv.template end<0>(); ++it)
219 { // loop through each element
220 int id = gv.indexSet().index(*it);
221 auto stress = myModel.getStress(id);
222 double h = myModel.FailureCriteria(stress);
223 if(h > max) max = h;
224 }
225 MPI_Allreduce(&max, &mech_fail, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
226 if(gv.comm().rank() == 0) std::cout << "Pure mechanical failure " << mech_fail << std::endl;
227 }*/
228 //Correction to mechanical load
229 double lambda = 1.0; //(1-therm_fail)/mech_fail;
230 // 3. Rerun with mechanical and thermal load and double check
231 // TODO should probably introduce a loop here in case this is still quite far off
232 double total_fail = 0;
233 int cnt; //count number of solves
234 while(std::abs(total_fail-1.0) > 1e-2){ //tolerance to which to determine failure load
235 if(gv.comm().rank() == 0) std::cout << "Lambda " << lambda << std::endl;
236 myModel.setThermal(true);
237 myModel.setPressure(lambda);
238
239 typedef Dune::PDELab::linearelasticity<GV, MODEL, dofel> LOP;
240 LOP lop(gv, myModel);
241
242 typedef Dune::PDELab::GridOperator<GFS,GFS,LOP,MBE,RF,RF,RF,C,C> GO;
243 GO go(gfs,cg,gfs,cg,lop,mbe);
244
245 // === Make coefficent vector and initialize it from a function
246 typedef Dune::PDELab::Backend::Vector<GFS,double> V;
247 V u(gfs,0.0);
248 Dune::PDELab::interpolate(initial_solution,gfs,u);
249
250 // === Set non constrained dofs to zero
251 Dune::PDELab::set_shifted_dofs(cg,0.0,u);
252
253 // === Solve the linear system
254 myModel.template solve<GO,V,GFS,C,Constraints,MBE,LOP>(go,u,gfs,cg,constraints,mbe,lop);
255 cnt++;
256 //if(gv.comm().rank() == 0) std::cout << "Check solution does something! ||x||_2 = " << u.two_norm() << std::endl;
257
258 if(gv.comm().rank() == 0) std::cout << "=== Calculate Stresses" << std::endl;
259 calculateStresses<MODEL,V,GV,GFS,MBE>(myModel,u,gv,gfs,mbe);
260
261 // Pure thermal failure
262 using Dune::PDELab::Backend::native;
263 double max = 0.0;
264 typedef typename GV::Traits::template Codim<0>::Iterator ElementIterator;
265 for (ElementIterator it = gv.template begin<0>(); it!=gv.template end<0>(); ++it)
266 { // loop through each element
267 int id = gv.indexSet().index(*it);
268 auto stress = myModel.getStress(id);
269 double h = myModel.FailureCriteria(stress);
270 if(h > max) max = h;
271 }
272 MPI_Allreduce(&max, &total_fail, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
273 if(gv.comm().rank() == 0) std::cout << "Total failure " << total_fail << std::endl;
274
275 lambda = lambda/total_fail;
276
277 if(std::abs(total_fail-1.0)<1e-2){ //If done, plot and postprocess
278 if(gv.comm().rank() == 0) std::cout << "Number of solves " << cnt << std::endl;
279
280 Dune::SubsamplingVTKWriter<GV> vtkwriter(gv,0);
281 Dune::PDELab::addSolutionToVTKWriter(vtkwriter,gfs,u);
282 vtkwriter.write(myModel.vtk_displacement_output,Dune::VTK::appendedraw);
283
284 plotProperties<MODEL,V,GV,GFS,MBE>(myModel,gv,gfs,mbe);
285
286 myModel.template postprocess<GO,V,GFS,C,MBE,GV>(go,u,gfs,cg,gv,mbe);
287 }
288 }
289
290 }
291 } // end apply
292
293 };
294 }
295}
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::heatStaticDriver
Static driver including a thermal loading.
Definition: heatStaticDriver.hh:19
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 (Sep 4, 22:38, 2025)