Loading [MathJax]/extensions/tex2jax.js

dune-composites (unstable)

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;
53 YGRID yaspgrid(myModel.coords,myModel.periodic,myModel.overlap,helper.getCommunicator(),&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 typedef Dune::PDELab::ConformingDirichletConstraints CON_EXT; //needed only for Geneo solver
108
109 times[2] = watch.elapsed();
110 watch.reset();
111
112 if(elem_order == 2){ // == Element Order
113 const int element_order = 2;
114 const int dofel = 3 * 20;
115 const int non_zeros = 81;
116
117 if(helper.rank() == 0)
118 std::cout << "Piecewise quadratic serendipity elements" << std::endl;
119
120 typedef Dune::PDELab::SerendipityLocalFiniteElementMap<GV,typename GV::Grid::ctype,double,element_order> FEM;
121 FEM fem(gv);
122
123 using SCALAR_VBE = Dune::PDELab::ISTL::VectorBackend<>;
124 typedef Dune::PDELab::GridFunctionSpace<GV, FEM, CON, SCALAR_VBE> SCALAR_GFS;
125 typedef Dune::PDELab::GridFunctionSpace<GV, FEM, CON_EXT, SCALAR_VBE> SCALAR_GFS_EXT;
126 SCALAR_GFS U1(gv,fem,con); U1.name("U1"); SCALAR_GFS_EXT U1_EXT(gv,fem);
127 SCALAR_GFS U2(gv,fem,con); U2.name("U2"); SCALAR_GFS_EXT U2_EXT(gv,fem);
128 SCALAR_GFS U3(gv,fem,con); U3.name("U3"); SCALAR_GFS_EXT U3_EXT(gv,fem);
129
130 // Note that Vectors are blocked by Dune::PDELab::EntityBlockedOrderingTag
131 const int block_size = 3;
132 using VBE = Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed>;
133 using Ordering = Dune::PDELab::EntityBlockedOrderingTag;
134 typedef Dune::PDELab::PowerGridFunctionSpace <SCALAR_GFS,block_size,VBE,Ordering> GFS;
135 const GFS gfs(U1,U2,U3);
136 typedef Dune::PDELab::PowerGridFunctionSpace <SCALAR_GFS_EXT,block_size,VBE,Ordering> GFS_EXT;
137 const GFS_EXT gfs_ext(U1_EXT,U2_EXT,U3_EXT);
138
139 typedef typename GFS::template ConstraintsContainer<RF>::Type C;
140
141 // Make constraints map and initialize it from a function
142 C cg;
143 cg.clear();
144
145 BC U1_cc(gv,myModel,0), U2_cc(gv,myModel,1), U3_cc(gv,myModel,2);
146 Constraints constraints(U1_cc,U2_cc,U3_cc);
147 Dune::PDELab::constraints(constraints,gfs,cg);
148
149 MBE mbe(non_zeros); // Maximal number of nonzeroes per row
150
151 //Correction to mechanical load
152 double lambda = 1.0; //(1-therm_fail)/mech_fail;
153 // 3. Rerun with mechanical and thermal load and double check
154 // TODO should probably introduce a loop here in case this is still quite far off
155 double total_fail = 0;
156 int cnt; //count number of solves
157 while(std::abs(total_fail-1.0) > 1e-2){ //tolerance to which to determine failure load
158 if(gv.comm().rank() == 0) std::cout << "Lambda " << lambda << std::endl;
159 myModel.setThermal(true);
160 myModel.setPressure(lambda);
161
162 typedef Dune::PDELab::linearelasticity<GV, MODEL, dofel> LOP;
163 LOP lop(gv, myModel);
164
165 typedef Dune::PDELab::GridOperator<GFS,GFS,LOP,MBE,RF,RF,RF,C,C> GO;
166 GO go(gfs,cg,gfs,cg,lop,mbe);
167 typedef Dune::PDELab::GridOperator<GFS_EXT,GFS_EXT,LOP,MBE,RF,RF,RF,C,C> GO_EXT;
168
169 // === Make coefficent vector and initialize it from a function
170 typedef Dune::PDELab::Backend::Vector<GFS,double> V;
171 V u(gfs,0.0);
172 Dune::PDELab::interpolate(initial_solution,gfs,u);
173
174 // === Set non constrained dofs to zero
175 Dune::PDELab::set_shifted_dofs(cg,0.0,u);
176
177 // === Solve the linear system
178 myModel.template solve<GO,GO_EXT,V,GFS,GFS_EXT,C,Constraints,MBE,LOP>(go,u,gfs,gfs_ext,cg,constraints,mbe,lop);
179 cnt++;
180 //if(gv.comm().rank() == 0) std::cout << "Check solution does something! ||x||_2 = " << u.two_norm() << std::endl;
181
182 if(gv.comm().rank() == 0) std::cout << "=== Calculate Stresses" << std::endl;
183 calculateStresses<MODEL,V,GV,GFS,MBE>(myModel,u,gv,gfs,mbe);
184
185 // Pure thermal failure
186 using Dune::PDELab::Backend::native;
187 double max = 0.0;
188 typedef typename GV::Traits::template Codim<0>::Iterator ElementIterator;
189 for (ElementIterator it = gv.template begin<0>(); it!=gv.template end<0>(); ++it)
190 { // loop through each element
191 int id = gv.indexSet().index(*it);
192 auto stress = myModel.getStress(id);
193 double h = myModel.FailureCriteria(stress);
194 if(h > max) max = h;
195 }
196 MPI_Allreduce(&max, &total_fail, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
197 if(gv.comm().rank() == 0) std::cout << "Total failure " << total_fail << std::endl;
198
199 lambda = lambda/total_fail;
200
201 if(std::abs(total_fail-1.0)<1e-2){ //If done, plot and postprocess
202 if(gv.comm().rank() == 0) std::cout << "Number of solves " << cnt << std::endl;
203
204 Dune::SubsamplingVTKWriter<GV> vtkwriter(gv,Dune::refinementLevels(0));
205 Dune::PDELab::addSolutionToVTKWriter(vtkwriter,gfs,u);
206 vtkwriter.write(myModel.vtk_displacement_output,Dune::VTK::appendedraw);
207
208 plotProperties<MODEL,V,GV,GFS,MBE>(myModel,gv,gfs,mbe);
209
210 myModel.template postprocess<GO,V,GFS,C,MBE,GV>(go,u,gfs,cg,gv,mbe);
211 }
212 }
213
214 }
215 } // end apply
216
217 };
218 }
219}
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 (Apr 3, 22:46, 2025)