dune-composites (unstable)

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