dune-composites (2.5.1)

baseStructuredGridModel.hh
1#ifndef CREATE_STRUCTUREMODEL__HH
2#define CREATE_STRUCTUREMODEL__HH
3
4#include <dune/composites/Setup/Cijkl.hh>
5#include <dune/composites/Setup/Region.hh>
6
7#include <dune/composites/Driver/Solvers/solver_info.hh>
8#include <dune/composites/Driver/Solvers/solvers.hh>
9
10using namespace std;
11
12namespace Dune{
13 namespace Composites{
14
16
20 class baseStructuredGridModel{
21
22 public:
23 //Types
24 typedef Dune::YaspGrid<3,Dune::TensorProductCoordinates<double,3>> YGRID;
25 typedef YGRID::LeafGridView YGV;
26 typedef Dune::FieldVector<double,3> FieldVec;
27 typedef Dune::FieldVector<double,6> Tensor2;
28 typedef Dune::FieldMatrix<double,6,6> Tensor4;
29
31 int overlap;
32
33 std::array<std::vector<double>,3> coords;
34
36 std::bitset<3> periodic;
37 Dune::MPIHelper& helper;
38 std::vector<Tensor4> baseTensor;
39
41 int numMaterials;
42 std::vector<Material> myMaterials;
43 std::vector<Region> R;
44
45 std::vector<int> elemIndx2PG;
46 std::vector<double> rot123;
47 std::vector<Tensor4> C;
48
49 bool verbosity;
50 double tolerance;
51 int solver_type;
52 int solver_tolerance;
53 int Krylov_Verb;
54 int Krylov_Max_It;
55 int solver_verb;
56
57 bool residualHeat;
58
59 std::vector<double> stackSeq;
60 std::vector<std::vector<double>> sig;
61 bool stress_Plot;
62 bool sigResized;
63 std::string vtk_displacement_output;
64 std::string vtk_stress_output;
65 std::string vtk_properties_output;
66 int intOrder;
67 std::vector<int> nel;
68 bool isParallel;
69 bool setUp_Required;
70
71 SolverInfo solverParameters;
72
73 double q; //pressure value
74
75 bool failureCriterionDefined;
76
77 int refineBaseGrid;
78
79 double gravity = 9.80665;
81
84 baseStructuredGridModel(Dune::MPIHelper& helper_, bool GenEO = false):helper(helper_){
85
86 isParallel = false;
87 if (helper.size() > 1){
88 isParallel = true;
89 }
90
91 /*solver_type = config.get<int>("solver.solver_type", 1);
92 Krylov_Max_It = 1000; //config.get<int>("solver.Krylov_Max_It", 500);
93 Krylov_Verb = config.get<int>("solver.Krylov_Verb", 0);
94 Krylov_Tol = config.get<double>("solver.Krylov_Tol", 1e-4);*/
95 intOrder = config.get<int>("fe.intOrder",5);
96
97 solver_verb = 1;
98 sigResized = false;
99 stress_Plot = true;
100
101 vtk_displacement_output = "baseModel_displacement";
102 vtk_stress_output = "baseModel_stress";
103 vtk_properties_output = "baseModel_properties";
104
105 residualHeat = false;
106 setUp_Required = true;
107 failureCriterionDefined = false;
108 refineBaseGrid = 0;
109 };
110
111 void inline setup(){
112 // setup -
113 LayerCake();
114 setUpMaterials();
115 computeElasticTensors();
116 }
117
118
119 void inline sizeStressVector(int nel){
120 sig.resize(nel);
121 for (unsigned int i = 0; i < sig.size(); i++){
122 sig[i].resize(6);
123 }
124 }
125
126 void inline setStress(int id, Dune::FieldVector<double,6>& element_stress){
127 for (int i = 0; i < 6; i++){
128 sig[id][i] = element_stress[i];
129 }
130 }
131
132 Dune::FieldVector<double,6> inline getStress(int id){
133 Dune::FieldVector<double,6> element_stress;
134 for (int i = 0; i < 6; i++){
135 element_stress[i] = sig[id][i];
136 }
137 return element_stress;
138 }
139
140 int inline getIntOrder() const {
141 return 5;
142 }
143
145 void inline setThermal(bool t) {
146 residualHeat = t;
147 }
148
149 int inline whichRegion(const FieldVec& x){
150 double Lx = 0.0;
151 int ans = 0;
152 bool flag = false;
153 for (unsigned int k = 0; k < R.size(); k++){
154 Lx += R[k].L[0];
155 if (x[0] < Lx && flag == false){ ans = k; flag = true; }
156 }
157 return ans;
158 }
159
160 int inline whichLayer(FieldVec& x, unsigned int r){
161 assert(r < R.size());
162 bool flag = false;
163 double z = 0.0;
164 int ans = 0;
165 for (unsigned int k = 0; k < R[r].layers.size(); k++){
166 z += R[r].layers[k].getThickness();
167 if (x[2] < z && flag == false){ ans = k; flag = true; }
168 }
169 return ans;
170 }
171
173 void inline evaluateHeat(Tensor2& f,int id) const{
174 f = 0.0;
175 if(residualHeat != false){
176 double deltaT = -160.; //temperature difference
177 double alpha11 = -0.342 * 1e-6;
178 double alpha22 = 25.8 * 1e-6;
179 double alpha33 = 25.8 * 1e-6;
180 double alpharesin = 25.8 * 1e-6;
181 if(getMaterialTypeFromElement(id) == 0){
182 f[0] = alpharesin*deltaT;
183 f[1] = alpharesin*deltaT;
184 f[2] = alpharesin*deltaT;
185 }
186 else{
187 f[0] = alpha11*deltaT;
188 f[1] = alpha22*deltaT;
189 f[2] = alpha33*deltaT;
190 }
191 }
192 }
193
194 inline void evaluateWeight(Dune::FieldVector<double,3>& f, int id) const{
195 int pg = elemIndx2PG[id];
196 f[2] = - gravity * myMaterials[pg].density;
197 }
198
199 inline void setPressure(double q_new) {
200 q = q_new;
201 }
202
204 inline void evaluateNeumann(const Dune::FieldVector<double,3> &x,Dune::FieldVector<double,3>& h,const Dune::FieldVector<double,3>& normal) const{
205 h = 0.0;
206 }
207
209 void inline setSolver(){
210 solver_type = 0; // UmfPack as default solver
211 if (helper.size() > 1){ // Parallel Solver
212 solver_type = 1;
213 }
214 }
215
216 int inline getElements(int i){
217 return nel[i];
218 }
219
228 void inline LayerCake(){
229
230 // Defines Default values for Base Class - can be overwritten by defining your own derived class
231
232 // Define Periodic Boundary Conditions
233 periodic[0] = false; periodic[1] = false; periodic[2] = false;
234
235 // Uniform Tensor Product Grid on [0,1]^3
236 std::vector<int> nel(3);
237 nel[0] = 5; nel[1] = 5; nel[2] = 5;
238 for (int i = 0; i < 3; i++){
239 double h = 1. / nel[i];
240 coords[i].resize(nel[i] + 1);
241 coords[i][0] = 0.0;
242 for (int j = 1; j <= nel[i]; j++){
243 coords[i][j] = coords[i][j-1] + h;
244 }
245 }
246 }
247
256 void inline LayerCakeFromFile(string pathToCsv){
257
258 ifstream csvFile(pathToCsv);
259 if(! csvFile) throw std::runtime_error("Could not open file " + pathToCsv);
260
261 string item;
262 int numRegions = 0;
263 int maxPly = 0;
264
265 string line;
266 getline(csvFile, line);
267 istringstream lineStream(line);
268
269 for (int j = 0; j < 2; j++){
270 // Obtain number of regions
271 getline(lineStream, item, ',');
272 if(helper.rank() == 0) std::cout << item << std::endl;
273 }
274
275 numRegions = atof(item.c_str());
276 if(helper.rank() == 0) std::cout << "Number of Regions = " << numRegions << std::endl;
277 R.resize(numRegions);
278
279 getline(csvFile, line);
280 istringstream lS2(line);
281
282 // Obtain max number of plies
283 getline(lS2, item, ',');
284 getline(lS2, item, ',');
285 maxPly = atof(item.c_str());
286 if(helper.rank() == 0) std::cout << "maxPly = " << maxPly << std::endl;
287
288 getline(csvFile, line);
289 istringstream lS3(line);
290 for (int j = 0; j < 4; j++){
291 // Obtain max number of plies
292 getline(lS3, item, ',');
293 if (j > 0){
294 periodic[j-1] = atof(item.c_str());
295 }
296 }
297 if(helper.rank() == 0) std::cout << "periodic = " << periodic << std::endl;
298
299 for(int j = 0; j < numRegions; j++){ // For each regions
300 //if(helper.rank() == 0) std::cout << "========= Region = " << j << std::endl;
301
302 R[j].setNumLayers(maxPly);
303 getline(csvFile, line); getline(csvFile, line);
304 int numParameters = 1;
305 istringstream lS4(line);
306 for (int ii = 0; ii < numParameters; ii++){
307 getline(lS4, item, ',');
308 if (ii == 0){
309 R[j].setType(atof(item.c_str()));
310 numParameters = R[j].numParameters;
311 }
312 else{
313 R[j].setParameter(ii-1,atof(item.c_str()));
314 }
315 } // For each parameter
316
317 for (int k = 0; k < maxPly; k++){
318 //if(helper.rank() == 0) std::cout << "Ply Number = " << k << std::endl;
319 getline(csvFile, line); // Obtain next row
320 istringstream lS5(line);
321 getline(lS5,item,',');
322 R[j].layers[k].setElements(atof(item.c_str()));
323
324 //if(helper.rank() == 0) std::cout << "Elements " << atof(item.c_str()) << std::endl;
325 getline(lS5,item,',');
326
327 //if(helper.rank() == 0) std::cout << "Material " << atof(item.c_str()) << std::endl;
328 R[j].layers[k].setMaterial(atof(item.c_str()) );
329 getline(lS5,item,',');
330 R[j].layers[k].setThickness(atof(item.c_str()));
331
332 //if(helper.rank() == 0) std::cout << "Thickness " << atof(item.c_str()) << std::endl;
333 getline(lS5,item,',');
334 R[j].layers[k].setOrientation(atof(item.c_str()));
335
336 //if(helper.rank() == 0) std::cout << "Orientation " << atof(item.c_str()) << std::endl;
337
338 } // for each ply
339 } // For each region
340 }
341
342 void inline GeometryBuilder(){
343 // Builds Overall Geometry and Mesh from a set of regions
344 if(helper.rank() == 0) std::cout << "=== Building Geometry" << std::endl;
345
346 Dune::FieldVector<double,3> origin(0.0); // Define Origin, the default of this will be set to [0.0,0.0,0.0]
347
348 // Compute number of elements in x and y. Since structures are current prismatic and tensor product grid y elements must be the same so we check this
349 int numRegions = R.size();
350
351 int nelx = 0;
352 int nely = R[0].nel[1];
353
354 for (int i = 0; i < numRegions;i++){
355 nelx += R[i].nel[0];
356 assert(R[i].nel[1] == nely); //"Overall tensor product grid y must be equal across all regions, check geometry input file");
357 }
358
359 if(helper.rank() == 0) std::cout << "Elements in x direction " << std::endl;
360
361 coords[0].resize(nelx + 1);
362 coords[0][0] = origin[0];
363
364 int k = 1; // initialise counter
365
366 for (int i = 0; i < numRegions; i++){
367 double hx = R[i].L[0] / R[i].nel[0]; // This makes grid uniform across a region TO DO:: Update for non-uniform grids
368 for (int j = 0; j < R[i].nel[0]; j++){
369 coords[0][k] = coords[0][k-1] + hx;
370 k++; // Increment Counter in the x-direction
371 }
372 }
373
374 // === Prismatic mesh in y
375 coords[1].resize(nely + 1);
376 coords[1][0] = origin[0];
377 double hy = R[0].L[1] / R[0].nel[1];
378 for (int j = 1; j < R[0].nel[1] + 1; j++){
379 coords[1][j] = coords[1][j-1] + hy;
380 }
381
382 // Build Coordinates in z
383 // For now let's assume uniform thickness (not necessary though)
384 int nelz = 0;
385 for (unsigned int i = 0; i < R[0].layers.size(); i++){
386 // for each layer
387 nelz += R[0].layers[i].getElements();
388 }
389 coords[2].resize(nelz + 1);
390 coords[2][0] = origin[2];
391
392 int e = 0;
393 for (unsigned int i = 0; i < R[0].layers.size(); i++){
394 double hz = R[0].layers[i].getThickness() / R[0].layers[i].getElements();
395 for (int k = 0; k < R[0].layers[i].getElements(); k++){
396 e++;
397 coords[2][e] = coords[2][e-1] + hz;
398 }
399 }
400 R[0].L[2] = coords[2][e];
401 }
402
403
404 int inline getSizeCoords(int i){ return coords[i].size(); }
405
406
407 void setUpMaterials(){
408 numMaterials = 3;
409 myMaterials.resize(numMaterials);
410
411 // Material I: Isotropic Resin
412 myMaterials[0].setType(0); // Isotropic
413 myMaterials[0].setProp(0,10000.); // - E, Young's Modulus
414 myMaterials[0].setProp(1,0.35); // - nu, Poisson Ratio
415 myMaterials[0].setDensity(1.57e-5); // - rho, Density kg/mm^2
416
417 // Material II: Orthotropic, Composite Ply - AS4/8552
418 myMaterials[1].setType(1); // Orthotropic Composite
419 myMaterials[1].setProp(0,162000.); // E11
420 myMaterials[1].setProp(1,10000.); // E22
421 myMaterials[1].setProp(2,10000.); // E33
422 myMaterials[1].setProp(3,0.35); // nu12
423 myMaterials[1].setProp(4,0.35); // nu13
424 myMaterials[1].setProp(5,0.49); // nu23
425 myMaterials[1].setProp(6,5200.); // G12
426 myMaterials[1].setProp(7,5200.); // G13
427 myMaterials[1].setProp(8,3500.); // G23
428 myMaterials[1].setDensity(1.57e-5); // - rho, Density kg/mm^2
429
430 // Material I: Isotropic Resin
431 myMaterials[2].setType(0); // Isotropic
432 myMaterials[2].setProp(0,10000000.); // - E, Young's Modulus
433 myMaterials[2].setProp(1,0.35); // - nu, Poisson Ratio
434 myMaterials[2].setDensity(1.57e-5); // - rho, Density kg/mm^2
435 }
436
437
438 double inline FailureCriteria(Dune::FieldVector<double,6>& stress) const{
439 return stress[0];
440 }
441
442 template <class GV, class GV2>
443 void computeTensorsOnGrid(const GV& gv, const GV2& gv2){
444
445 typedef typename GV::Traits::template Codim<0>::Iterator ElementIterator;
446 auto is = gv.indexSet();
447 C.resize(gv.size(0));
448
449 // Calculate Cijkl[id] in each element
450 auto it2 = gv2.template begin<0>();
451
452 for (ElementIterator it = gv.template begin<0>(); it!=gv.template end<0>(); ++it)
453 { // loop through each element
454 int id = is.index(*it); // Get element id
455 int pg = elemIndx2PG[id];
456 C[id] = TensorRotation(baseTensor[pg],rot123[id],2);
457 it2++;
458 } // End for-loop of each element */
459 }; // End Constructor
460
461 FieldVec inline gridTransformation(const FieldVec & x) const{
462 // Default grid transformation is the identity map
463 auto y = x;
464 return y;
465 }
466
467 std::array<int, 3> inline GridPartition(){
468 // Function returns default partition
469 std::array<int, 3> Partitioning;
470 std::fill(Partitioning.begin(),Partitioning.end(), 1);
471 return Partitioning;
472 }
473
474 void inline setPG(const YGV& ygv){
475 elemIndx2PG.resize(ygv.size(0));
476 rot123.resize(ygv.size(0));
477 // Loop Through Each element and Assign Physical Group and Rotation
478 for (const auto& eit : elements(ygv)){
479 int id = ygv.indexSet().index(eit);
480 auto point = eit.geometry().center();
481 //elemIndx2PG[id] = setPhysicalGroup(point);
482 // Which Layer?
483 int theRegion = whichRegion(point);
484 int theLayer = whichLayer(point,theRegion);
485 elemIndx2PG[id] = R[theRegion].layers[theLayer].getMaterial();
486 rot123[id] = R[theRegion].layers[theLayer].getOrientation();
487 for(int i = 0; i < eit.geometry().corners(); i++){ //Loop over nodes of the element
488 auto pc = eit.geometry().corner(i);
489 if(isMPC(pc)){ //Check for mpc constraints
490 elemIndx2PG[id] = 2;
491 rot123[id] = 0;
492 }
493 }
494 }
495 }
496
497 int inline getMaterialType(int whichMaterial){
498 assert(whichMaterial < numMaterials);
499 return myMaterials[whichMaterial].Type;
500 }
501
502
503 int inline setPhysicalGroup(FieldVec& x){
504 return 0;
505 }
506
507 void inline computeElasticTensors(){
508 // Calculate Elastic Tensors for each material
509 baseTensor.resize(numMaterials);
510 for (int i = 0; i < numMaterials; i++){
511 baseTensor[i] = Cijkl(myMaterials[i]);
512 }
513 }
514
515
516 Tensor4 TensorRotation(Tensor4& baseC, double theta, int i_dim) const
517 {
518 Tensor4 locC = baseC;
519 Tensor4 R2(0.0),R2T(0.0);
520
521 //Set up 3D rotation matrix
522 double c = cos(theta*M_PI/180.0);
523 double s = sin(theta*M_PI/180.0);
524 Dune::FieldMatrix<double,3,3> A;
525 if(i_dim == 2) A = {{c,s,0},{-s,c,0},{0,0,1}};
526 if(i_dim == 1) A = {{c,0,-s},{0,1,0},{s,0,c}};
527 if(i_dim == 0) A = {{1,0,0},{0,c,s},{0,-s,c}};
528
529 //Transform to tensor rotation and transpose
530 R2 = getMatrix(A);
531 transpose<Tensor4,6>(R2,R2T);
532
533 //rotate the input matrix
534 locC.leftmultiply(R2);
535 locC.rightmultiply(R2T);
536 return locC;
537 } // end eval
538
539 Tensor4 TensorRotationLeft(Tensor4& C_local, double theta, int i_dim) const{
540 Tensor4 locC = C_local;
541 Tensor4 R2(0.0),R2T(0.0);
542
543 //Set up 3D rotation matrix
544 double c = cos(theta*M_PI/180.0);
545 double s = sin(theta*M_PI/180.0);
546 Dune::FieldMatrix<double,3,3> A;
547 if(i_dim == 2) A = {{c,s,0},{-s,c,0},{0,0,1}};
548 if(i_dim == 1) A = {{c,0,-s},{0,1,0},{s,0,c}};
549 if(i_dim == 0) A = {{1,0,0},{0,c,s},{0,-s,c}};
550
551 //Transform to tensor rotation
552 R2 = getMatrix(A);
553
554 //rotate the input matrix
555 locC.leftmultiply(R2);
556 return locC;
557 }
558
559 Tensor4 TensorRotationRight(Tensor4& C_local, double theta, int i_dim) const{
560 Tensor4 locC = C_local;
561 Tensor4 R2(0.0),R2T(0.0);
562
563 //set up 3D rotation matrix
564 double c = cos(theta*M_PI/180.0);
565 double s = sin(theta*M_PI/180.0);
566 Dune::FieldMatrix<double,3,3> A;
567 if(i_dim == 2) A = {{c,s,0},{-s,c,0},{0,0,1}};
568 if(i_dim == 1) A = {{c,0,-s},{0,1,0},{s,0,c}};
569 if(i_dim == 0) A = {{1,0,0},{0,c,s},{0,-s,c}};
570
571 //Transform to tensor rotation and transpose
572 R2 = getMatrix(A);
573 transpose<Tensor4,6>(R2,R2T);
574
575 //rotate the input matrix
576 locC.rightmultiply(R2T);
577 return locC;
578 }
579
580 private:
581 //Turn 3D rotation matrix into a tensor rotation
582 Tensor4 getMatrix(Dune::FieldMatrix<double,3,3> A) const{
583 Tensor4 R;
584 R[0][0] = A[0][0]*A[0][0]; R[0][1] = A[0][1]*A[0][1]; R[0][2] = A[0][2]*A[0][2];
585 R[1][0] = A[1][0]*A[1][0]; R[1][1] = A[1][1]*A[1][1]; R[1][2] = A[1][2]*A[1][2];
586 R[2][0] = A[2][0]*A[2][0]; R[2][1] = A[2][1]*A[2][1]; R[2][2] = A[2][2]*A[2][2];
587
588 R[0][3] = 2.0*A[0][1]*A[0][2]; R[0][4] = 2.0*A[0][0]*A[0][2]; R[0][5] = 2.0*A[0][0]*A[0][1];
589 R[1][3] = 2.0*A[1][1]*A[1][2]; R[1][4] = 2.0*A[1][0]*A[1][2]; R[1][5] = 2.0*A[1][0]*A[1][1];
590 R[2][3] = 2.0*A[2][1]*A[2][2]; R[2][4] = 2.0*A[2][0]*A[2][2]; R[2][5] = 2.0*A[2][0]*A[2][1];
591
592 R[3][0] = A[1][0]*A[2][0]; R[3][1] = A[1][1]*A[2][1]; R[3][2] = A[1][2]*A[2][2];
593 R[4][0] = A[0][0]*A[2][0]; R[4][1] = A[0][1]*A[2][1]; R[4][2] = A[0][2]*A[2][2];
594 R[5][0] = A[0][0]*A[1][0]; R[5][1] = A[0][1]*A[1][1]; R[5][2] = A[0][2]*A[1][2];
595
596 R[3][3] = A[1][1]*A[2][2]+A[1][2]*A[2][1]; R[3][4] = A[1][0]*A[2][2]+A[1][2]*A[2][0]; R[3][5] = A[1][0]*A[2][1]+A[1][1]*A[2][0];
597 R[4][3] = A[0][1]*A[2][2]+A[2][1]*A[0][2]; R[4][4] = A[0][0]*A[2][2]+A[2][0]*A[0][2]; R[4][5] = A[2][0]*A[0][1]+A[2][1]*A[0][0];
598 R[5][3] = A[0][1]*A[1][2]+A[0][2]*A[1][1]; R[5][4] = A[0][0]*A[1][2]+A[0][2]*A[1][0]; R[5][5] = A[0][0]*A[1][1]+A[0][1]*A[1][0];
599 return R;
600 }
601
602 //void inline evaluateTensor(int id, Tensor4& C_){ C_ = C[id]; }
603 public:
604
605 Tensor4 inline getElasticTensor(int id) const { return C[id]; }
606
607 std::vector<double> returnTheta(const FieldVec& x) const{
608 FieldVec y, yph;
609 double h = 1e-10;
610 FieldVec xph = x; xph[0] += h/2.;
611 FieldVec xmh = x; xmh[0] -= h/2.;
612 y = gridTransformation(xph);
613 yph = gridTransformation(xmh);
614 auto angle = atan2((yph[2]-y[2])/h,(yph[0]-y[0])/h);
615 std::vector<double> theta(3);
616 theta[1] = angle;
617 return theta;
618 }
619
622 bool inline isDirichlet(FieldVec& x, int dof){
623 // Default isDirichlet function, returning homogeneous Dirichlet boundary conditioners
624 if (x[0] < 1e-6)
625 return true;
626 return false;
627 }
628
630 bool inline isMPC(FieldVec& x){
631 return false;
632 }
633
635 double inline evaluateDirichlet(const FieldVec& x, int dof) const{
636 // Return homogeneous boundary conditions
637 return 0.0;
638 }
639
640 double inline getMaterialTypeFromElement(int id) const{
641 int pg = elemIndx2PG[id];
642 return myMaterials[pg].Type;
643 }
644
645 double inline getOrientation(int id) const {return rot123[id];}
646
647 template<class GO, class U, class GFS, class C,class MBE, class GV>
648 void inline postprocess(const GO& go, U& u, const GFS& gfs, const C& cg, const GV gv, MBE& mbe){
649 // Default does nothing
650 }
651
652 template<class GO, class V, class GFS, class C, class CON, class MBE, class LOP>
653 bool inline solve(GO& go, V& u, const GFS& gfs, C& cg, CON& constraints, MBE& mbe, LOP& lop){
654
655 bool flag = true;
656
657 if(gfs.gridView().comm().size() > 1){ // == DEFAULT AVAILABLE PARALLEL SOLVERS
658 if(solverParameters.preconditioner == "GenEO"){
659 typedef Dune::Geneo::CG_GenEO<V,GFS,GO,LOP,CON,MBE> GenEO;
660 int cells = std::pow(gfs.gridView().size(0),1./3.);
661 GenEO mysolver(u,gfs,lop,constraints,mbe,solverParameters,helper, cells);
662 mysolver.apply();
663 return flag;
664 }
665 if(solverParameters.preconditioner == "OneLevel_AdditiveSchwarz"){
666 typedef Dune::PDELab::ISTLBackend_OVLP_CG_UMFPack<GFS, C> PAR_UMF; //CG with UMFPack
667 PAR_UMF ls(gfs,cg,solverParameters.MaxIt,solverParameters.verb);
668 Dune::PDELab::StationaryLinearProblemSolver<GO,PAR_UMF,V> slp(go,ls,u,solverParameters.KrylovTol);
669 slp.apply();
670 return flag;
671 }
672 }
673 if(solverParameters.solver == "Hypre"){
674#if HAVE_HYPRE
675 typedef typename GO::Jacobian Matrix;
676 Matrix a(go);
677 auto b = u;
678 go.residual(u,b);
679 b *= -1.0;
680 HypreSolver<GFS, Matrix, V> solver(gfs,u,a,b, solverParameters.KrylovTol);
681 go.jacobian(u,a);
682 solver.solve(u);
683 return flag;
684#endif
685 }
686 if(gfs.gridView().comm().size() == 1){ // == SEQUENTIAL SOLVER
687 if(solverParameters.solver == "UMFPack"){
688#if HAVE_SUITESPARSE
689 typedef Dune::PDELab::ISTLBackend_SEQ_UMFPack SEQ_UMF;
690 SEQ_UMF ls(solver_verb); //Direct solver
691 Dune::PDELab::StationaryLinearProblemSolver<GO, SEQ_UMF,V> slp(go,ls,u,1e-6); // note tolerance is ignored as Direct solver
692 slp.apply();
693#else
694 if(gfs.gridView().rank() == 0) std::cout << "SEQ_UMFPack solver not available" << std::endl;
695#endif
696 return flag;
697 }
698 if(solverParameters.preconditioner == "AMG"){
699 typedef Dune::PDELab::ISTLBackend_SEQ_CG_AMG_SSOR<GO> SEQ_CG_AMG_SSOR;
700 SEQ_CG_AMG_SSOR ls(500,verbosity);
701 Dune::PDELab::StationaryLinearProblemSolver<GO,SEQ_CG_AMG_SSOR,V> slp(go,ls,u,tolerance);
702 slp.apply();
703 }
704 }
705
706 return flag;
707 }
708
709
710 private:
711 bool GenEO;
712
713 };
714
715 }
716}
717
718#endif
Dune::Composites::baseStructuredGridModel
structuredGridModel class.
Definition: baseStructuredGridModel.hh:20
Dune::Composites::baseStructuredGridModel::solver_tolerance
int solver_tolerance
Variable which controls solver type (UMFPack for sequential and GENEO for parallel)
Definition: baseStructuredGridModel.hh:52
Dune::Composites::baseStructuredGridModel::isMPC
bool isMPC(FieldVec &x)
A function which returns whether a point is constrained by a multi-point-constraint.
Definition: baseStructuredGridModel.hh:630
Dune::Composites::baseStructuredGridModel::LayerCake
void LayerCake()
A member taking which constructs the base layered laminate in untransformed coordinates.
Definition: baseStructuredGridModel.hh:228
Dune::Composites::baseStructuredGridModel::Krylov_Verb
int Krylov_Verb
Tolerance of Solver (ignored for a direct solver)
Definition: baseStructuredGridModel.hh:53
Dune::Composites::baseStructuredGridModel::rot123
std::vector< double > rot123
std::vector containing information of each material defined
Definition: baseStructuredGridModel.hh:46
Dune::Composites::baseStructuredGridModel::verbosity
bool verbosity
Vector storing elastic tensor within each element.
Definition: baseStructuredGridModel.hh:49
Dune::Composites::baseStructuredGridModel::setThermal
void setThermal(bool t)
returns true when thermal strain is applied
Definition: baseStructuredGridModel.hh:145
Dune::Composites::baseStructuredGridModel::Krylov_Max_It
int Krylov_Max_It
Variable Controlling Verbosity of Krylov Solver.
Definition: baseStructuredGridModel.hh:54
Dune::Composites::baseStructuredGridModel::LayerCakeFromFile
void LayerCakeFromFile(string pathToCsv)
Function to read LayerCake from File.
Definition: baseStructuredGridModel.hh:256
Dune::Composites::baseStructuredGridModel::R
std::vector< Region > R
std::vector containing information of each material defined
Definition: baseStructuredGridModel.hh:43
Dune::Composites::baseStructuredGridModel::nel
std::vector< int > nel
Gaussian order of integration (default value is 5)
Definition: baseStructuredGridModel.hh:67
Dune::Composites::baseStructuredGridModel::periodic
std::bitset< 3 > periodic
Array which stores periodicity of grid in each of the 3 dimensions, default value are set to all fals...
Definition: baseStructuredGridModel.hh:36
Dune::Composites::baseStructuredGridModel::solver_type
int solver_type
Tolerance of iterative solver.
Definition: baseStructuredGridModel.hh:51
Dune::Composites::baseStructuredGridModel::evaluateNeumann
void evaluateNeumann(const Dune::FieldVector< double, 3 > &x, Dune::FieldVector< double, 3 > &h, const Dune::FieldVector< double, 3 > &normal) const
function which evaluates the Neumann boundary conditions at the point x
Definition: baseStructuredGridModel.hh:204
Dune::Composites::baseStructuredGridModel::q
double q
Instance of SolverInfo to define solver parameters and store solver results e.g. iterations,...
Definition: baseStructuredGridModel.hh:73
Dune::Composites::baseStructuredGridModel::overlap
int overlap
Size of overlap between subdomains for domain decomposition.
Definition: baseStructuredGridModel.hh:31
Dune::Composites::baseStructuredGridModel::setUp_Required
bool setUp_Required
Boolean variable which indicates if run is parallel or not, this later will allow to run multiple sol...
Definition: baseStructuredGridModel.hh:69
Dune::Composites::baseStructuredGridModel::elemIndx2PG
std::vector< int > elemIndx2PG
std::vector containing information of each Region
Definition: baseStructuredGridModel.hh:45
Dune::Composites::baseStructuredGridModel::solver_verb
int solver_verb
Variable define max iterations.
Definition: baseStructuredGridModel.hh:55
Dune::Composites::baseStructuredGridModel::numMaterials
int numMaterials
Number of Different types of materials defined for applications.
Definition: baseStructuredGridModel.hh:41
Dune::Composites::baseStructuredGridModel::C
std::vector< Tensor4 > C
Vector defining rotation of ply in material coordinates.
Definition: baseStructuredGridModel.hh:47
Dune::Composites::baseStructuredGridModel::sigResized
bool sigResized
Variable if stress is plotted out in vtk format.
Definition: baseStructuredGridModel.hh:62
Dune::Composites::baseStructuredGridModel::evaluateHeat
void evaluateHeat(Tensor2 &f, int id) const
Function which evaluates the thermal stress tensor sig_thermal on an element.
Definition: baseStructuredGridModel.hh:173
Dune::Composites::baseStructuredGridModel::evaluateDirichlet
double evaluateDirichlet(const FieldVec &x, int dof) const
returns the Dirichlet boundary value at a point for the degree of freedom dof.
Definition: baseStructuredGridModel.hh:635
Dune::Composites::baseStructuredGridModel::isParallel
bool isParallel
Number of elements in each direction.
Definition: baseStructuredGridModel.hh:68
Dune::Composites::baseStructuredGridModel::tolerance
double tolerance
Variable controls level of output to interfacea.
Definition: baseStructuredGridModel.hh:50
Dune::Composites::baseStructuredGridModel::baseStructuredGridModel
baseStructuredGridModel(Dune::MPIHelper &helper_, bool GenEO=false)
The main constructor.
Definition: baseStructuredGridModel.hh:84
Dune::Composites::baseStructuredGridModel::isDirichlet
bool isDirichlet(FieldVec &x, int dof)
Definition: baseStructuredGridModel.hh:622
Dune::Composites::baseStructuredGridModel::stackSeq
std::vector< double > stackSeq
Set whether to use residual heat stresses, default false.
Definition: baseStructuredGridModel.hh:59
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden & Uni Heidelberg  |  generated with Hugo v0.111.3 (Sep 5, 22:35, 2025)