Finite Element Domain Decomposition Library
FEDDLib
Loading...
Searching...
No Matches
FE_ElementAssembly_def.hpp
1#ifndef FE_ELEMENTASSEMBLY_DEF_hpp
2#define FE_ELEMENTASSEMBLY_DEF_hpp
3
4#include <string>
5#include "feddlib/core/core_config.h"
6#ifdef FEDD_HAVE_ACEGENINTERFACE
7#include <aceinterface.hpp>
8#endif
9
18
19
20
21using Teuchos::reduceAll;
22using Teuchos::REDUCE_SUM;
23using Teuchos::outArg;
24
25namespace FEDD {
26
27
28template <class SC, class LO, class GO, class NO>
29FE_ElementAssembly<SC,LO,GO,NO>::FE_ElementAssembly(bool saveAssembly):
30domainVec_(0),
31setZeros_(false),
32myeps_(),
33saveAssembly_(saveAssembly)
34{
35}
36
37template <class SC, class LO, class GO, class NO>
38void FE_ElementAssembly<SC,LO,GO,NO>::addFE(DomainConstPtr_Type domain){
39
40 if (saveAssembly_){
41 DomainPtr_Type domainNC = Teuchos::rcp_const_cast<Domain_Type>( domain );
42 domainNC->initializeFEData();
43 }
44 domainVec_.push_back(domain);
45
46}
47
48template <class SC, class LO, class GO, class NO>
49void FE_ElementAssembly<SC,LO,GO,NO>::doSetZeros(double eps){
50
51 setZeros_ = true;
52 myeps_ = eps;
53
54}
55
67
68// Check the order of chemistry and solid in system matrix
69/*template <class SC, class LO, class GO, class NO>
70void FE_ElementAssembly<SC,LO,GO,NO>::globalAssembly(std::string ProblemType,
71 int dim,
72 int degree,
73 MultiVectorPtr_Type sol_rep,
74 BlockMatrixPtr_Type &A,
75 BlockMultiVectorPtr_Type &resVec,
76 ParameterListPtr_Type params,
77 std::string assembleMode,
78 bool callFillComplete,
79 int FELocExternal){
80
81 int problemSize = domainVec_.size(); // Problem size, as in number of subproblems
82
83 // Depending on problem size we extract necessary information
84
87 int numChem=3;
88 if(FETypeChem == "P2"){
89 numChem=6;
90 }
91 if(dim==3){
92 numChem=4;
93 if(FETypeChem == "P2")
94 numChem=10;
95 }
96 int numSolid=3;
97 if(FETypeSolid == "P2")
98 numSolid=6;
99
100 if(dim==3){
101 numSolid=4;
102 if(FETypeSolid == "P2")
103 numSolid=10;
104 }
105
106 BlockMultiVectorPtr_Type resVecRep = Teuchos::rcp( new BlockMultiVector_Type( 2) );
107
108 tuple_disk_vec_ptr_Type problemDisk = Teuchos::rcp(new tuple_disk_vec_Type(0));
109 for(int i=0; i< problemSize; i++){
110 int numNodes=0.;
111 if(domainVec.at(i)->getFEType() == "P2"){
112 numNodes=6;
113 }
114 if(dim==3){
115 numNodes=4;
116 if(domainVec.at(i)->getFEType() == "P2")
117 numNodes=10;
118 }
119 tuple_ssii_Type infoTuple (domainVec.at(i)->getPhysic(),domainVec.at(i)->getFEType(),domainVec.at(i)->getDofs(),numChem);
120 problemDisk->push_back(infoTuple);
121
122 MultiVectorPtr_Type resVec;
123 if(domainVec.at(i)->getDofs() == 1)
124 resVec = Teuchos::rcp( new MultiVector_Type( domainVec_.at(i)->getMapRepeated(), 1 ) );
125 else
126 resVec = Teuchos::rcp( new MultiVector_Type( domainVec_.at(i)->getMapVecFieldRepeated(), 1 ) );
127
128 resVecRep->addBlock(resVec,i);
129
130
131 }
132
133
134 if(assemblyFEElements_.size()== 0){
135 initAssembleFEElements(ProblemType,problemDisk,domainVec_.at(0)->getElementsC(), params,domainVec_.at(0)->getPointsRepeated(),domainVec_.at(0)->getElementMap());
136 }
137 else if(assemblyFEElements_.size() != domainVec_.at(0)->getElementsC()->numberElements())
138 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Number Elements not the same as number assembleFE elements." );
139
140 vec_dbl_Type solution_tmp;
141 for (UN T=0; T<assemblyFEElements_.size(); T++) {
142 vec_dbl_Type solution(0);
143 vec_FE_Type elements;
144 for(int i=0; i< problemSize; i++){
145 solution_tmp = getSolution(domainVec_.at(i)->getElementsC()->getElement(T).getVectorNodeList(), sol_rep->getBlock(i),domainVec.at(i)->getDofs());
146 solution.insert( solution.end(), solution_tmp.begin(), solution_tmp.end() );
147
148
149 }
150
151 assemblyFEElements_[T]->updateSolution(solution);
152
153 SmallMatrixPtr_Type elementMatrix;
154 vec_dbl_ptr_Type rhsVec;
155
156 if(assembleMode == "Jacobian"){
157 assemblyFEElements_[T]->assembleJacobian();
158
159 elementMatrix = assemblyFEElements_[T]->getJacobian();
160 //elementMatrix->print();
161 assemblyFEElements_[T]->advanceNewtonStep(); // n genereal non linear solver step
162
163 addFeBlockMatrix(A, elementMatrix, domainVec_.at(0)->getElementsC()->getElement(T), problemDisk);
164
165 }
166 if(assembleMode == "Rhs"){
167 assemblyFEElements_[T]->assembleRHS();
168 rhsVec = assemblyFEElements_[T]->getRHS();
169 addFeBlockMv(resVecRep, rhsVec, elementsSolid->getElement(T),elementsChem->getElement(T), dofsSolid,dofsChem);
170
171 }
172
173 }
174 if ( assembleMode == "Jacobian"){
175 for(int probRow = 0; probRow < problemSize; probRow++){
176 for(int probCol = 0; probCol < problemSize; probCol++){
177 MapConstPtr_Type rowMap;
178 MapConstPtr_Type domainMap;
179
180 if(domainVec.at(probRow)->getDofs() == 1)
181 rowMap = domainVec_.at(FElocRow)->getMapUnique();
182 else
183 rowMap = domainVec_.at(probRow)->getMapVecFieldUnique();
184
185 if(domainVec.at(probCol)->getDofs() == 1)
186 domainMap = domainVec_.at(probCol)->getMapUnique()
187 else
188 domainMap = domainVec_.at(probCol)->getMapVecFieldUnique()
189
190 A->getBlock(probRow,probCol)->fillComplete(domainMap,rowMap);
191 }
192 }
193 }
194 else if(assembleMode == "Rhs"){
195
196 for(int probRow = 0; probRow < problemSize; probRow++){
197 MapConstPtr_Type rowMap;
198 if(domainVec.at(probRow)->getDofs() == 1)
199 rowMap = domainVec_.at(FElocRow)->getMapUnique();
200 else
201 rowMap = domainVec_.at(probRow)->getMapVecFieldUnique();
202
203 MultiVectorPtr_Type resVecUnique = Teuchos::rcp( new MultiVector_Type( rowMap, 1 ) );
204
205 resVecUnique->putScalar(0.);
206
207 resVecUnique->exportFromVector( resVecRep, true, "Add" );
208
209 resVec->addBlock(resVecUnique,probRow);
210 }
211 }
212
213}*/
214
217 \brief Inserting element stiffness matrices into global stiffness matrix
218
219
220@param[in] &A Global Block Matrix
221@param[in] elementMatrix Stiffness matrix of one element
222@param[in] element Corresponding finite element
223@param[in] map Map that corresponds to repeated nodes of first block
224@param[in] map Map that corresponds to repeated nodes of second block
225
226*/
227/*template <class SC, class LO, class GO, class NO>
228void FE_ElementAssembly<SC,LO,GO,NO>::addFeBlockMatrix(BlockMatrixPtr_Type &A, SmallMatrixPtr_Type elementMatrix, FiniteElement_vec element, tuple_disk_vec_ptr_Type problemDisk){
229
230 int numDisk = problemDisk->size();
231
232 for(int probRow = 0; probRow < numDisk; probRow++){
233 for(int probCol = 0; probCol < numDisk; probCol++){
234 int dofs1 = std::get<2>(problemDisk->at(probRow));
235 int dofs2 = std::get<2>(problemDisk->at(probCol));
236
237 int numNodes1 = std::get<3>(problemDisk->at(probRow));
238 int numNodes2=std::get<3>(problemDisk->at(probCol));
239
240 int offsetRow= numNodes1*dofs1*probRow;
241 int offsetCol= numNodes1*dofs1*probCol;
242
243 mapRow = domainVec.at(probRow)->getMapRepeated();
244 mapCol = domainVec.at(probCol)->getMapRepeated();
245
246 Teuchos::Array<SC> value2( numNodes2, 0. );
247 Teuchos::Array<GO> columnIndices2( numNodes2, 0 );
248 for (UN i=0; i < numNodes1; i++){
249 for(int di=0; di<dofs1; di++){
250 GO row =GO (dofs1* mapRow->getGlobalElement( element[probRow].getNode(i) )+di);
251 for(int d=0; d<dofs2; d++){
252 for (UN j=0; j < numNodes2 ; j++) {
253 value2[j] = (*elementMatrix)[offsetRow+i*dofs1+di][offsetCol+j*dofs2+d];
254 columnIndices2[j] =GO (dofs2* mapCol->getGlobalElement( element.getNode(j) )+d);
255 }
256 A->getBlock(probRow,probCol)->insertGlobalValues( row, columnIndices2(), value2() ); // Automatically adds entries if a value already exists
257 }
258 }
259 }
260
261 }
262 }
263}
264
265*/
266
277/*
278template <class SC, class LO, class GO, class NO>
279void FE_ElementAssembly<SC,LO,GO,NO>::addFeBlockMv(BlockMultiVectorPtr_Type &res, vec_dbl_ptr_Type rhsVec,tuple_disk_vec_ptr_Type problemDisk, FiniteElement_vec element)
280 int numDisk = problemDisk->size();
281
282 for(int probRow = 0; probRow < numDisk; probRow++){
283 Teuchos::ArrayRCP<SC> resArray_block = res->getBlockNonConst(probRow)->getDataNonConst(0);
284 vec_LO_Type nodeList_block = element[probRow].getVectorNodeList();
285 int dofs = std::get<2>(problemDisk->at(probRow));
286 int numNodes = std::get<3>(problemDisk->at(probRow));
287
288 int offset = numNodes*dofs;
289 for(int i=0; i< nodeList_block.size() ; i++){
290 for(int d=0; d<dofs; d++){
291 resArray_block[nodeList_block[i]*dofs+d] += (*rhsVec)[i*dofs+d+offset];
292 }
293 }
294 }
295}
296*/
297
309
310template <class SC, class LO, class GO, class NO>
311void FE_ElementAssembly<SC,LO,GO,NO>::assemblyLinearElasticity(int dim,
312 std::string FEType,
313 int degree,
314 int dofs,
315 MultiVectorPtr_Type d_rep,
316 BlockMatrixPtr_Type &A,
317 BlockMultiVectorPtr_Type &resVec,
318 ParameterListPtr_Type params,
319 bool reAssemble,
320 std::string assembleMode,
321 bool callFillComplete,
322 int FELocExternal){
323
324 ElementsPtr_Type elements = domainVec_.at(0)->getElementsC();
325
326 int dofsElement = elements->getElement(0).getVectorNodeList().size();
327
328 vec2D_dbl_ptr_Type pointsRep = domainVec_.at(0)->getPointsRepeated();
329
330 MapConstPtr_Type mapVel = domainVec_.at(0)->getMapRepeated();
331
332 vec_dbl_Type solution(0);
333 vec_dbl_Type solution_d;
334
335 vec_dbl_Type rhsVec;
336
339 int numNodes=6;
340 if(dim==3){
341 numNodes=10;
342 }
343 tuple_disk_vec_ptr_Type problemDisk = Teuchos::rcp(new tuple_disk_vec_Type(0));
344 tuple_ssii_Type displacement ("Displacement",FEType,dofs,numNodes);
345 problemDisk->push_back(displacement);
346
347 if(assemblyFEElements_.size()== 0)
348 initAssembleFEElements("LinearElasticity",problemDisk,elements, params,pointsRep,domainVec_.at(0)->getElementMap());
349 else if(assemblyFEElements_.size() != elements->numberElements())
350 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Number Elements not the same as number assembleFE elements." );
351
352 MultiVectorPtr_Type resVec_d = Teuchos::rcp( new MultiVector_Type( domainVec_.at(0)->getMapVecFieldRepeated(), 1 ) );
353
354 BlockMultiVectorPtr_Type resVecRep = Teuchos::rcp( new BlockMultiVector_Type( 1) );
355 resVecRep->addBlock(resVec_d,0);
356
357 SmallMatrixPtr_Type elementMatrix;
358 for (UN T=0; T<assemblyFEElements_.size(); T++) {
359 vec_dbl_Type solution(0);
360
361 solution_d = getSolution(elements->getElement(T).getVectorNodeList(), d_rep,dofs);
362
363 solution.insert( solution.end(), solution_d.begin(), solution_d.end() );
364
365 assemblyFEElements_[T]->updateSolution(solution);
366
367 assemblyFEElements_[T]->assembleJacobian();
368
369 elementMatrix = assemblyFEElements_[T]->getJacobian();
370
371 assemblyFEElements_[T]->advanceNewtonStep();
372
373
374 addFeBlock(A, elementMatrix, elements->getElement(T), mapVel, 0, 0, problemDisk);
375
376 }
377 if (callFillComplete)
378 A->getBlock(0,0)->fillComplete( domainVec_.at(0)->getMapVecFieldUnique(),domainVec_.at(0)->getMapVecFieldUnique());
379
380
381
382
383}
384
396
397template <class SC, class LO, class GO, class NO>
398void FE_ElementAssembly<SC,LO,GO,NO>::assemblyNonLinearElasticity(int dim,
399 std::string FEType,
400 int degree,
401 int dofs,
402 MultiVectorPtr_Type d_rep,
403 BlockMatrixPtr_Type &A,
404 BlockMultiVectorPtr_Type &resVec,
405 ParameterListPtr_Type params,
406 bool callFillComplete,
407 int FELocExternal){
408
409 ElementsPtr_Type elements = domainVec_.at(0)->getElementsC();
410
411 int dofsElement = elements->getElement(0).getVectorNodeList().size();
412
413 vec2D_dbl_ptr_Type pointsRep = domainVec_.at(0)->getPointsRepeated();
414
415 MapConstPtr_Type map = domainVec_.at(0)->getMapRepeated();
416
417 vec_dbl_Type solution(0);
418 vec_dbl_Type solution_d;
419
420 vec_dbl_ptr_Type rhsVec;
421
424 int numNodes=6;
425 if(dim==3){
426 numNodes=10;
427 }
428 tuple_disk_vec_ptr_Type problemDisk = Teuchos::rcp(new tuple_disk_vec_Type(0));
429 tuple_ssii_Type displacement ("Displacement",FEType,dofs,numNodes);
430 problemDisk->push_back(displacement);
431
432 int neoHookeNum = params->sublist("Parameter").get("Neo-Hooke Modell",1);
433
434 std::string nonLinElasModell = "NonLinearElasticity2";
435 if(neoHookeNum == 1)
436 nonLinElasModell = "NonLinearElasticity";
437
438 //std::cout << " ######## Assembly Modell: " << nonLinElasModell << " ############ " << std::endl;
439
440
441 if(assemblyFEElements_.size()== 0)
442 initAssembleFEElements(nonLinElasModell,problemDisk,elements, params,pointsRep,domainVec_.at(0)->getElementMap());
443 else if(assemblyFEElements_.size() != elements->numberElements())
444 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Number Elements not the same as number assembleFE elements." );
445
446
447 SmallMatrixPtr_Type elementMatrix;
448 for (UN T=0; T<assemblyFEElements_.size(); T++) {
449 vec_dbl_Type solution(0);
450
451 solution_d = getSolution(elements->getElement(T).getVectorNodeList(), d_rep,dofs);
452
453 solution.insert( solution.end(), solution_d.begin(), solution_d.end() );
454
455 assemblyFEElements_[T]->updateSolution(solution);
456
457 assemblyFEElements_[T]->assembleJacobian();
458 elementMatrix = assemblyFEElements_[T]->getJacobian();
459 addFeBlock(A, elementMatrix, elements->getElement(T), map, 0, 0, problemDisk);
460
461 assemblyFEElements_[T]->assembleRHS();
462 rhsVec = assemblyFEElements_[T]->getRHS();
463 addFeBlockMv(resVec, rhsVec, elements->getElement(T), dofs);
464
465 assemblyFEElements_[T]->advanceNewtonStep();
466
467
468 }
469 if (callFillComplete)
470 A->getBlock(0,0)->fillComplete( domainVec_.at(0)->getMapVecFieldUnique(),domainVec_.at(0)->getMapVecFieldUnique());
471
472}
473
485
486template <class SC, class LO, class GO, class NO>
487void FE_ElementAssembly<SC,LO,GO,NO>::assemblyNonLinearElasticity(int dim,
488 std::string FEType,
489 int degree,
490 int dofs,
491 MultiVectorPtr_Type d_rep,
492 BlockMatrixPtr_Type &A,
493 BlockMultiVectorPtr_Type &resVec,
494 ParameterListPtr_Type params,
495 DomainConstPtr_Type domain,
496 MultiVectorPtr_Type eModVec,
497 bool callFillComplete,
498 int FELocExternal){
499
500 ElementsPtr_Type elements = domain->getElementsC();
501
502 int dofsElement = elements->getElement(0).getVectorNodeList().size();
503
504 vec2D_dbl_ptr_Type pointsRep = domain->getPointsRepeated();
505
506 MapConstPtr_Type map = domain->getMapRepeated();
507
508 vec_dbl_Type solution(0);
509 vec_dbl_Type solution_d;
510
511 vec_dbl_ptr_Type rhsVec;
512
515 int numNodes=6;
516 if(dim==3){
517 numNodes=10;
518 }
519 tuple_disk_vec_ptr_Type problemDisk = Teuchos::rcp(new tuple_disk_vec_Type(0));
520 tuple_ssii_Type displacement ("Displacement",FEType,dofs,numNodes);
521 problemDisk->push_back(displacement);
522
523 int neoHookeNum = params->sublist("Parameter").get("Neo-Hooke Modell",1);
524
525 std::string nonLinElasModell = "NonLinearElasticity2";
526 if(neoHookeNum == 1)
527 nonLinElasModell = "NonLinearElasticity";
528
529 //std::cout << " ######## Assembly Modell: " << nonLinElasModell << " ############ " << std::endl;
530
531 if(assemblyFEElements_.size()== 0)
532 initAssembleFEElements(nonLinElasModell,problemDisk,elements, params,pointsRep,domain->getElementMap());
533 else if(assemblyFEElements_.size() != elements->numberElements())
534 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Number Elements not the same as number assembleFE elements." );
535
536 Teuchos::ArrayRCP<SC> eModVecA = eModVec->getDataNonConst(0);
537
538 SmallMatrixPtr_Type elementMatrix;
539 for (UN T=0; T<assemblyFEElements_.size(); T++) {
540 vec_dbl_Type solution(0);
541
542 solution_d = getSolution(elements->getElement(T).getVectorNodeList(), d_rep,dofs);
543
544 solution.insert( solution.end(), solution_d.begin(), solution_d.end() );
545
546 assemblyFEElements_[T]->updateSolution(solution);
547 assemblyFEElements_[T]->updateParameter("E",eModVecA[T]);
548 assemblyFEElements_[T]->assembleJacobian();
549 elementMatrix = assemblyFEElements_[T]->getJacobian();
550 addFeBlock(A, elementMatrix, elements->getElement(T), map, 0, 0, problemDisk);
551
552 assemblyFEElements_[T]->assembleRHS();
553 rhsVec = assemblyFEElements_[T]->getRHS();
554 addFeBlockMv(resVec, rhsVec, elements->getElement(T), dofs);
555
556 assemblyFEElements_[T]->advanceNewtonStep();
557
558
559 }
560 if (callFillComplete)
561 A->getBlock(0,0)->fillComplete( domainVec_.at(0)->getMapVecFieldUnique(),domainVec_.at(0)->getMapVecFieldUnique());
562
563}
564
565
576
577template <class SC, class LO, class GO, class NO>
578void FE_ElementAssembly<SC,LO,GO,NO>::addFeBlockMv(BlockMultiVectorPtr_Type &res, vec_dbl_ptr_Type rhsVec, FiniteElement elementBlock, int dofs){
579
580 Teuchos::ArrayRCP<SC> resArray_block = res->getBlockNonConst(0)->getDataNonConst(0);
581
582 vec_LO_Type nodeList_block = elementBlock.getVectorNodeList();
583
584 for(int i=0; i< nodeList_block.size() ; i++){
585 for(int d=0; d<dofs; d++)
586 resArray_block[nodeList_block[i]*dofs+d] += (*rhsVec)[i*dofs+d];
587 }
588}
589
590// Check the order of chemistry and solid in system matrix
591template <class SC, class LO, class GO, class NO>
592void FE_ElementAssembly<SC,LO,GO,NO>::assemblyAceDeformDiffu(int dim,
593 std::string FETypeChem,
594 std::string FETypeSolid,
595 int degree,
596 int dofsChem,
597 int dofsSolid,
598 MultiVectorPtr_Type c_rep,
599 MultiVectorPtr_Type d_rep,
600 BlockMatrixPtr_Type &A,
601 BlockMultiVectorPtr_Type &resVec,
602 ParameterListPtr_Type params,
603 std::string assembleMode,
604 bool callFillComplete,
605 int FELocExternal){
606
607 if((FETypeChem != "P2") || (FETypeSolid != "P2") || dim != 3)
608 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "No AceGen Implementation available for Discretization and Dimension." );
609
610
611 UN FElocChem = 1; //checkFE(dim,FETypeChem); // Checks for different domains which belongs to a certain fetype
612 UN FElocSolid = 0; //checkFE(dim,FETypeSolid); // Checks for different domains which belongs to a certain fetype
613
614 ElementsPtr_Type elementsChem= domainVec_.at(FElocChem)->getElementsC();
615
616 ElementsPtr_Type elementsSolid = domainVec_.at(FElocSolid)->getElementsC();
617
618 //this->domainVec_.at(FElocChem)->info();
619 //this->domainVec_.at(FElocSolid)->info();
620 //int dofsElement = elements->getElement(0).getVectorNodeList().size();
621
622 vec2D_dbl_ptr_Type pointsRep = domainVec_.at(FElocSolid)->getPointsRepeated();
623
624 MapConstPtr_Type mapChem = domainVec_.at(FElocChem)->getMapRepeated();
625
626 MapConstPtr_Type mapSolid = domainVec_.at(FElocSolid)->getMapRepeated();
627
628 vec_dbl_Type solution_c;
629 vec_dbl_Type solution_d;
630
631 vec_dbl_ptr_Type rhsVec;
632
635 int numChem=3;
636 if(FETypeChem == "P2"){
637 numChem=6;
638 }
639 if(dim==3){
640 numChem=4;
641 if(FETypeChem == "P2")
642 numChem=10;
643 }
644 int numSolid=3;
645 if(FETypeSolid == "P2")
646 numSolid=6;
647
648 if(dim==3){
649 numSolid=4;
650 if(FETypeSolid == "P2")
651 numSolid=10;
652 }
653 tuple_disk_vec_ptr_Type problemDisk = Teuchos::rcp(new tuple_disk_vec_Type(0));
654 tuple_ssii_Type chem ("Chemistry",FETypeChem,dofsChem,numChem);
655 tuple_ssii_Type solid ("Solid",FETypeSolid,dofsSolid,numSolid);
656 problemDisk->push_back(solid);
657 problemDisk->push_back(chem);
658
659 tuple_disk_vec_ptr_Type problemDiskChem = Teuchos::rcp(new tuple_disk_vec_Type(0));
660 problemDiskChem->push_back(chem);
661
662 std::string SCIModel = params->sublist("Parameter").get("Structure Model","SCI_simple");
663
664 if(assemblyFEElements_.size()== 0){
665 initAssembleFEElements(SCIModel,problemDisk,elementsChem, params,pointsRep,domainVec_.at(FElocSolid)->getElementMap());
666 }
667 else if(assemblyFEElements_.size() != elementsChem->numberElements())
668 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Number Elements not the same as number assembleFE elements." );
669
670 //SmallMatrixPtr_Type elementMatrix =Teuchos::rcp( new SmallMatrix_Type( dofsElement));
671
672 MultiVectorPtr_Type resVec_c = Teuchos::rcp( new MultiVector_Type( domainVec_.at(FElocChem)->getMapRepeated(), 1 ) );
673 MultiVectorPtr_Type resVec_d = Teuchos::rcp( new MultiVector_Type( domainVec_.at(FElocSolid)->getMapVecFieldRepeated(), 1 ) );
674
675 BlockMultiVectorPtr_Type resVecRep = Teuchos::rcp( new BlockMultiVector_Type( 2) );
676 resVecRep->addBlock(resVec_d,0);
677 resVecRep->addBlock(resVec_c,1);
678
679 SC detB;
680 SC absDetB;
681 SmallMatrix<SC> B(dim);
682 SmallMatrix<SC> Binv(dim);
683
684
685 for (UN T=0; T<assemblyFEElements_.size(); T++) {
686 vec_dbl_Type solution(0);
687
688 solution_c = getSolution(elementsChem->getElement(T).getVectorNodeList(), c_rep,dofsChem);
689 solution_d = getSolution(elementsSolid->getElement(T).getVectorNodeList(), d_rep,dofsSolid);
690
691 // First Solid, then Chemistry
692 solution.insert( solution.end(), solution_d.begin(), solution_d.end() );
693 solution.insert( solution.end(), solution_c.begin(), solution_c.end() );
694
695 assemblyFEElements_[T]->updateSolution(solution);
696
697 SmallMatrixPtr_Type elementMatrix;
698
699 // ------------------------
700 /*buildTransformation(elementsSolid->getElement(T).getVectorNodeList(), pointsRep, B, FETypeSolid);
701 detB = B.computeInverse(Binv);
702 absDetB = std::fabs(detB);
703 std::cout << " Determinante " << detB << std::endl;*/
704 // ------------------------
705
706
707
708
709 if(assembleMode == "Jacobian"){
710 assemblyFEElements_[T]->assembleJacobian();
711
712 elementMatrix = assemblyFEElements_[T]->getJacobian();
713 // elementMatrix->print();
714 assemblyFEElements_[T]->advanceNewtonStep(); // n genereal non linear solver step
715
716 addFeBlockMatrix(A, elementMatrix, elementsSolid->getElement(T), elementsSolid->getElement(T), mapSolid, mapChem, problemDisk);
717
718
719
720 }
721 if(assembleMode == "Rhs"){
722 assemblyFEElements_[T]->assembleRHS();
723 rhsVec = assemblyFEElements_[T]->getRHS();
724 addFeBlockMv(resVecRep, rhsVec, elementsSolid->getElement(T),elementsChem->getElement(T), dofsSolid,dofsChem);
725
726 }
727 if(assembleMode=="MassMatrix"){
728 assemblyFEElements_[T]->assembleJacobian();
729
730 AssembleFE_SCI_SMC_Active_Growth_Reorientation_Ptr_Type elTmp = Teuchos::rcp_dynamic_cast<AssembleFE_SCI_SMC_Active_Growth_Reorientation_Type>(assemblyFEElements_[T] );
731 elTmp->getMassMatrix(elementMatrix);
732 //elementMatrix->print();
733 addFeBlock(A, elementMatrix, elementsChem->getElement(T), mapChem, 0, 0, problemDiskChem);
734
735
736 }
737
738
739
740
741 }
742 if ( assembleMode == "Jacobian"){
743 A->getBlock(0,0)->fillComplete();
744 A->getBlock(1,0)->fillComplete(domainVec_.at(FElocSolid)->getMapVecFieldUnique(),domainVec_.at(FElocChem)->getMapUnique());
745 A->getBlock(0,1)->fillComplete(domainVec_.at(FElocChem)->getMapUnique(),domainVec_.at(FElocSolid)->getMapVecFieldUnique());
746 A->getBlock(1,1)->fillComplete();
747 }
748 else if(assembleMode == "Rhs"){
749
750 MultiVectorPtr_Type resVecUnique_d = Teuchos::rcp( new MultiVector_Type( domainVec_.at(FElocSolid)->getMapVecFieldUnique(), 1 ) );
751 MultiVectorPtr_Type resVecUnique_c = Teuchos::rcp( new MultiVector_Type( domainVec_.at(FElocChem)->getMapUnique(), 1 ) );
752
753 resVecUnique_d->putScalar(0.);
754 resVecUnique_c->putScalar(0.);
755
756 resVecUnique_d->exportFromVector( resVec_d, true, "Add" );
757 resVecUnique_c->exportFromVector( resVec_c, true, "Add" );
758
759 resVec->addBlock(resVecUnique_d,0);
760 resVec->addBlock(resVecUnique_c,1);
761 }
762 else if(assembleMode == "MassMatrix"){
763 A->getBlock(0,0)->fillComplete();
764 }
765
766}
767
768// Check the order of chemistry and solid in system matrix
769template <class SC, class LO, class GO, class NO>
770void FE_ElementAssembly<SC,LO,GO,NO>::assemblyAceDeformDiffuBlock(int dim,
771 std::string FETypeChem,
772 std::string FETypeSolid,
773 int degree,
774 int dofsChem,
775 int dofsSolid,
776 MultiVectorPtr_Type c_rep,
777 MultiVectorPtr_Type d_rep,
778 BlockMatrixPtr_Type &A,
779 int blockRow,
780 int blockCol,
781 BlockMultiVectorPtr_Type &resVec,
782 int block,
783 ParameterListPtr_Type params,
784 std::string assembleMode,
785 bool callFillComplete,
786 int FELocExternal){
787
788 if((FETypeChem != "P2") || (FETypeSolid != "P2") || dim != 3)
789 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "No AceGen Implementation available for Discretization and Dimension." );
790 if((blockRow != blockCol) && blockRow != 0)
791 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Block assemblyDeformDiffu AceGEN Version only implemented for 0,0 block right now" );
792
793 //UN FElocChem = 1; //checkFE(dim,FETypeChem); // Checks for different domains which belongs to a certain fetype
794 UN FElocSolid = 0; //checkFE(dim,FETypeSolid); // Checks for different domains which belongs to a certain fetype
795
796 ElementsPtr_Type elementsChem= domainVec_.at(FElocSolid)->getElementsC();
797
798 ElementsPtr_Type elementsSolid = domainVec_.at(FElocSolid)->getElementsC();
799
800 //this->domainVec_.at(FElocChem)->info();
801 //this->domainVec_.at(FElocSolid)->info();
802 //int dofsElement = elements->getElement(0).getVectorNodeList().size();
803
804 vec2D_dbl_ptr_Type pointsRep = domainVec_.at(FElocSolid)->getPointsRepeated();
805
806 //MapConstPtr_Type mapChem = domainVec_.at(FElocChem)->getMapRepeated();
807
808 MapConstPtr_Type mapSolid = domainVec_.at(FElocSolid)->getMapRepeated();
809
810 vec_dbl_Type solution_c;
811 vec_dbl_Type solution_d;
812
813 vec_dbl_ptr_Type rhsVec;
814
817 int numChem=3;
818 if(FETypeChem == "P2"){
819 numChem=6;
820 }
821 if(dim==3){
822 numChem=4;
823 if(FETypeChem == "P2")
824 numChem=10;
825 }
826 int numSolid=3;
827 if(FETypeSolid == "P2")
828 numSolid=6;
829
830 if(dim==3){
831 numSolid=4;
832 if(FETypeSolid == "P2")
833 numSolid=10;
834 }
835 tuple_disk_vec_ptr_Type problemDisk = Teuchos::rcp(new tuple_disk_vec_Type(0));
836 tuple_ssii_Type chem ("Chemistry",FETypeChem,dofsChem,numChem);
837 tuple_ssii_Type solid ("Solid",FETypeSolid,dofsSolid,numSolid);
838 problemDisk->push_back(solid);
839 problemDisk->push_back(chem);
840
841 std::string SCIModel = params->sublist("Parameter").get("Structure Model","SCI_simple");
842
843 if(assemblyFEElements_.size()== 0){
844 initAssembleFEElements(SCIModel,problemDisk,elementsChem, params,pointsRep,domainVec_.at(FElocSolid)->getElementMap());
845 }
846 else if(assemblyFEElements_.size() != elementsChem->numberElements())
847 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Number Elements not the same as number assembleFE elements." );
848
849 //SmallMatrixPtr_Type elementMatrix =Teuchos::rcp( new SmallMatrix_Type( dofsElement));
850
851 //MultiVectorPtr_Type resVec_c = Teuchos::rcp( new MultiVector_Type( domainVec_.at(FElocChem)->getMapRepeated(), 1 ) );
852 MultiVectorPtr_Type resVec_d = Teuchos::rcp( new MultiVector_Type( domainVec_.at(FElocSolid)->getMapVecFieldRepeated(), 1 ) );
853
854 BlockMultiVectorPtr_Type resVecRep = Teuchos::rcp( new BlockMultiVector_Type( 1) );
855 resVecRep->addBlock(resVec_d,0);
856 //resVecRep->addBlock(resVec_c,1);
857
858 for (UN T=0; T<assemblyFEElements_.size(); T++) {
859 vec_dbl_Type solution(0);
860
861 solution_c = getSolution(elementsChem->getElement(T).getVectorNodeList(), c_rep,dofsChem);
862 solution_d = getSolution(elementsSolid->getElement(T).getVectorNodeList(), d_rep,dofsSolid);
863
864 // First Solid, then Chemistry
865 solution.insert( solution.end(), solution_d.begin(), solution_d.end() );
866 solution.insert( solution.end(), solution_c.begin(), solution_c.end() );
867
868 assemblyFEElements_[T]->updateSolution(solution);
869
870 SmallMatrixPtr_Type elementMatrix;
871
872 if(assembleMode == "Jacobian"){
873 assemblyFEElements_[T]->assembleJacobian();
874
875 elementMatrix = assemblyFEElements_[T]->getJacobian();
876 // elementMatrix->print();
877 assemblyFEElements_[T]->advanceNewtonStep(); // n genereal non linear solver step
878 addFeBlock(A, elementMatrix, elementsSolid->getElement(T), mapSolid, 0, 0, problemDisk);
879 //addFeBlockMatrix(A, elementMatrix, elementsSolid->getElement(T), mapSolid, mapChem, problemDisk);
880 }
881 if(assembleMode == "Rhs"){
882 assemblyFEElements_[T]->assembleRHS();
883 rhsVec = assemblyFEElements_[T]->getRHS();
884 addFeBlockMv(resVecRep, rhsVec, elementsSolid->getElement(T), dofsSolid);
885 }
886 //if(assembleMode=="compute")
887 // assemblyFEElements_[T]->compute();
888
889
890
891
892 }
893 if ( assembleMode != "Rhs"){
894 A->getBlock(0,0)->fillComplete();
895 //A->getBlock(1,0)->fillComplete(domainVec_.at(FElocSolid)->getMapVecFieldUnique(),domainVec_.at(FElocChem)->getMapUnique());
896 //A->getBlock(0,1)->fillComplete(domainVec_.at(FElocChem)->getMapUnique(),domainVec_.at(FElocSolid)->getMapVecFieldUnique());
897 //A->getBlock(1,1)->fillComplete();
898 }
899
900 if(assembleMode == "Rhs"){
901
902 MultiVectorPtr_Type resVecUnique_d = Teuchos::rcp( new MultiVector_Type( domainVec_.at(FElocSolid)->getMapVecFieldUnique(), 1 ) );
903 //MultiVectorPtr_Type resVecUnique_c = Teuchos::rcp( new MultiVector_Type( domainVec_.at(FElocChem)->getMapUnique(), 1 ) );
904
905 resVecUnique_d->putScalar(0.);
906 //resVecUnique_c->putScalar(0.);
907
908 resVecUnique_d->exportFromVector( resVec_d, true, "Add" );
909 //resVecUnique_c->exportFromVector( resVec_c, true, "Add" );
910
911 resVec->addBlock(resVecUnique_d,0);
912 //resVec->addBlock(resVecUnique_c,1);
913 }
914
915
916}
917
929template <class SC, class LO, class GO, class NO>
930void FE_ElementAssembly<SC,LO,GO,NO>::addFeBlockMatrix(BlockMatrixPtr_Type &A, SmallMatrixPtr_Type elementMatrix, FiniteElement element1, FiniteElement element2, MapConstPtr_Type mapFirstRow,MapConstPtr_Type mapSecondRow, tuple_disk_vec_ptr_Type problemDisk){
931
932 int numDisk = problemDisk->size();
933
934 int dofs1 = std::get<2>(problemDisk->at(0));
935 int dofs2 = std::get<2>(problemDisk->at(1));
936
937 int numNodes1 = std::get<3>(problemDisk->at(0));
938 int numNodes2=std::get<3>(problemDisk->at(1));
939
940 int dofsBlock1 = dofs1*numNodes1;
941 int dofsBlock2 = dofs2*numNodes2;
942
943 Teuchos::Array<SC> value1( numNodes1, 0. );
944 Teuchos::Array<GO> columnIndices1( numNodes1, 0 );
945
946 Teuchos::Array<SC> value2( numNodes2, 0. );
947 Teuchos::Array<GO> columnIndices2( numNodes2, 0 );
948
949 for (UN i=0; i < numNodes1 ; i++) {
950 for(int di=0; di<dofs1; di++){
951 GO row =GO (dofs1* mapFirstRow->getGlobalElement( element1.getNode(i) )+di);
952 for(int d=0; d<dofs1; d++){
953 for (UN j=0; j < columnIndices1.size(); j++){
954 columnIndices1[j] = GO ( dofs1 * mapFirstRow->getGlobalElement( element1.getNode(j) ) + d );
955 value1[j] = (*elementMatrix)[dofs1*i+di][dofs1*j+d];
956 }
957 A->getBlock(0,0)->insertGlobalValues( row, columnIndices1(), value1() ); // Automatically adds entries if a value already exists
958 }
959 }
960 }
961 int offset= numNodes1*dofs1;
962
963
964 Teuchos::Array<SC> value( 1, 0. );
965 Teuchos::Array<GO> columnIndex( 1, 0 );
966 for (UN i=0; i < numNodes2 ; i++) {
967 for(int di=0; di<dofs2; di++){
968 GO row =GO (dofs2* mapSecondRow->getGlobalElement( element2.getNode(i) )+di);
969 for(int d=0; d<dofs2; d++){
970 for (UN j=0; j < columnIndices2.size(); j++){
971 double tmpValue = (*elementMatrix)[offset+dofs2*i+di][offset+dofs2*j+d];
972 if(std::fabs(tmpValue) > 1.e-13){
973 columnIndex[0] = GO ( dofs2 * mapSecondRow->getGlobalElement( element2.getNode(j) ) + d );
974 value[0] = tmpValue;
975 A->getBlock(1,1)->insertGlobalValues( row, columnIndex(), value() ); // Automatically adds entries if a value already exists
976
977 }
978 }
979 }
980 }
981 }
982
983 for (UN i=0; i < numNodes1; i++){
984 for(int di=0; di<dofs1; di++){
985 GO row =GO (dofs1* mapFirstRow->getGlobalElement( element1.getNode(i) )+di);
986 for(int d=0; d<dofs2; d++){
987 for (UN j=0; j < numNodes2 ; j++) {
988 value2[j] = (*elementMatrix)[i*dofs1+di][offset+j*dofs2+d];
989 columnIndices2[j] =GO (dofs2* mapSecondRow->getGlobalElement( element2.getNode(j) )+d);
990 }
991 A->getBlock(0,1)->insertGlobalValues( row, columnIndices2(), value2() ); // Automatically adds entries if a value already exists
992 }
993 }
994 }
995
996 for (UN j=0; j < numNodes2; j++){
997 for(int di=0; di<dofs2; di++){
998 GO row = GO (dofs2* mapSecondRow->getGlobalElement( element2.getNode(j) ) +di );
999 for(int d=0; d<dofs1; d++){
1000 for (UN i=0; i < numNodes1 ; i++) {
1001 value1[i] = (*elementMatrix)[offset+j*dofs2+di][dofs1*i+d];
1002 columnIndices1[i] =GO (dofs1* mapFirstRow->getGlobalElement( element1.getNode(i) )+d);
1003 }
1004 A->getBlock(1,0)->insertGlobalValues( row, columnIndices1(), value1() ); // Automatically adds entries if a value already exists
1005 }
1006 }
1007 }
1008
1009
1010}
1011
1012
1024
1025template <class SC, class LO, class GO, class NO>
1026void FE_ElementAssembly<SC,LO,GO,NO>::assemblyNavierStokes(int dim,
1027 std::string FETypeVelocity,
1028 std::string FETypePressure,
1029 int degree,
1030 int dofsVelocity,
1031 int dofsPressure,
1032 MultiVectorPtr_Type u_rep,
1033 MultiVectorPtr_Type p_rep,
1034 BlockMatrixPtr_Type &A,
1035 BlockMultiVectorPtr_Type &resVec,
1036 SmallMatrix_Type coeff,
1037 ParameterListPtr_Type params,
1038 bool reAssemble,
1039 std::string assembleMode,
1040 bool callFillComplete,
1041 int FELocExternal){
1042
1043
1044 UN FElocVel = checkFE(dim,FETypeVelocity); // Checks for different domains which belongs to a certain fetype
1045 UN FElocPres = checkFE(dim,FETypePressure); // Checks for different domains which belongs to a certain fetype
1046
1047 ElementsPtr_Type elements = domainVec_.at(FElocVel)->getElementsC();
1048
1049 ElementsPtr_Type elementsPres = domainVec_.at(FElocPres)->getElementsC();
1050
1051 int dofsElement = elements->getElement(0).getVectorNodeList().size();
1052
1053 vec2D_dbl_ptr_Type pointsRep = domainVec_.at(FElocVel)->getPointsRepeated();
1054
1055 MapConstPtr_Type mapVel = domainVec_.at(FElocVel)->getMapRepeated();
1056
1057 MapConstPtr_Type mapPres = domainVec_.at(FElocPres)->getMapRepeated();
1058
1059 vec_dbl_Type solution(0);
1060 vec_dbl_Type solution_u;
1061 vec_dbl_Type solution_p;
1062
1063 vec_dbl_ptr_Type rhsVec;
1064
1067 int numVelo=3;
1068 if(FETypeVelocity == "P2")
1069 numVelo=6;
1070
1071 if(dim==3){
1072 numVelo=4;
1073 if(FETypeVelocity == "P2")
1074 numVelo=10;
1075 }
1076 tuple_disk_vec_ptr_Type problemDisk = Teuchos::rcp(new tuple_disk_vec_Type(0));
1077 tuple_ssii_Type vel ("Velocity",FETypeVelocity,dofsVelocity,numVelo);
1078 tuple_ssii_Type pres ("Pressure",FETypePressure,dofsPressure,dim+1);
1079 problemDisk->push_back(vel);
1080 problemDisk->push_back(pres);
1081
1082 if(assemblyFEElements_.size()== 0){
1083 if(params->sublist("Material").get("Newtonian",true) == false)
1084 initAssembleFEElements("GeneralizedNewtonian",problemDisk,elements, params,pointsRep,domainVec_.at(FElocVel)->getElementMap()); // In cas of non Newtonian Fluid
1085 else
1086 initAssembleFEElements("NavierStokes",problemDisk,elements, params,pointsRep,domainVec_.at(FElocVel)->getElementMap());
1087 }
1088 else if(assemblyFEElements_.size() != elements->numberElements())
1089 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Number Elements not the same as number assembleFE elements." );
1090
1091 //SmallMatrixPtr_Type elementMatrix =Teuchos::rcp( new SmallMatrix_Type( dofsElement));
1092
1093 MultiVectorPtr_Type resVec_u = Teuchos::rcp( new MultiVector_Type( domainVec_.at(FElocVel)->getMapVecFieldRepeated(), 1 ) );
1094 MultiVectorPtr_Type resVec_p = Teuchos::rcp( new MultiVector_Type( domainVec_.at(FElocPres)->getMapRepeated(), 1 ) );
1095
1096 BlockMultiVectorPtr_Type resVecRep = Teuchos::rcp( new BlockMultiVector_Type( 2) );
1097 resVecRep->addBlock(resVec_u,0);
1098 resVecRep->addBlock(resVec_p,1);
1099
1100
1101 for (UN T=0; T<assemblyFEElements_.size(); T++) {
1102 vec_dbl_Type solution(0);
1103
1104 solution_u = getSolution(elements->getElement(T).getVectorNodeList(), u_rep,dofsVelocity);
1105 solution_p = getSolution(elementsPres->getElement(T).getVectorNodeList(), p_rep,dofsPressure);
1106
1107 solution.insert( solution.end(), solution_u.begin(), solution_u.end() );
1108 solution.insert( solution.end(), solution_p.begin(), solution_p.end() );
1109
1110 assemblyFEElements_[T]->updateSolution(solution);
1111
1112 SmallMatrixPtr_Type elementMatrix;
1113
1114 if(assembleMode == "Jacobian"){
1115 assemblyFEElements_[T]->assembleJacobian();
1116
1117 elementMatrix = assemblyFEElements_[T]->getJacobian();
1118
1119 assemblyFEElements_[T]->advanceNewtonStep(); // n genereal non linear solver step
1120
1121 if(reAssemble)
1122 addFeBlock(A, elementMatrix, elements->getElement(T), mapVel, 0, 0, problemDisk);
1123 else
1124 addFeBlockMatrix(A, elementMatrix, elements->getElement(T), elementsPres->getElement(T), mapVel, mapPres, problemDisk);
1125 }
1126 if(assembleMode == "FixedPoint"){
1127
1128 //AssembleFENavierStokesPtr_Type elTmp = Teuchos::rcp_dynamic_cast<AssembleFENavierStokes_Type>(assemblyFEElements_[T] ); // Why should we need a pointer we can directly call it using assemblyFEElements_[T]? Because assemblyFixedPoint is not a function of base class
1129
1130 if(params->sublist("Material").get("Newtonian",true) == false)
1131 {
1132 AssembleFEGeneralizedNewtonianPtr_Type elTmp = Teuchos::rcp_dynamic_cast<AssembleFEGeneralizedNewtonian_Type>( assemblyFEElements_[T] );
1133 elTmp->assembleFixedPoint();
1134 elementMatrix = elTmp->getFixedPointMatrix();
1135 }
1136 else // Newtonian Case
1137 {
1138 AssembleFENavierStokesPtr_Type elTmp = Teuchos::rcp_dynamic_cast<AssembleFENavierStokes_Type>( assemblyFEElements_[T] );
1139 elTmp->assembleFixedPoint();
1140 elementMatrix = elTmp->getFixedPointMatrix();
1141 }
1142
1143
1144 assemblyFEElements_[T]->advanceNewtonStep(); // n genereal non linear solver step
1145
1146 if(reAssemble)
1147 addFeBlock(A, elementMatrix, elements->getElement(T), mapVel, 0, 0, problemDisk);
1148 else
1149 addFeBlockMatrix(A, elementMatrix, elements->getElement(T), elementsPres->getElement(T),mapVel, mapPres, problemDisk);
1150
1151 }
1152 if(assembleMode == "Rhs"){
1153 AssembleFENavierStokesPtr_Type elTmp = Teuchos::rcp_dynamic_cast<AssembleFENavierStokes_Type>(assemblyFEElements_[T] );
1154 elTmp->setCoeff(coeff);// Coeffs from time discretization. Right now default [1][1] // [0][0]
1155 assemblyFEElements_[T]->assembleRHS();
1156 rhsVec = assemblyFEElements_[T]->getRHS();
1157 addFeBlockMv(resVecRep, rhsVec, elements->getElement(T),elementsPres->getElement(T), dofsVelocity,dofsPressure);
1158 }
1159
1160
1161 }
1162 if (callFillComplete && reAssemble && assembleMode != "Rhs" )
1163 A->getBlock(0,0)->fillComplete( domainVec_.at(FElocVel)->getMapVecFieldUnique(),domainVec_.at(FElocVel)->getMapVecFieldUnique());
1164 else if(callFillComplete && !reAssemble && assembleMode != "Rhs"){
1165 A->getBlock(0,0)->fillComplete();
1166 A->getBlock(1,0)->fillComplete(domainVec_.at(FElocVel)->getMapVecFieldUnique(),domainVec_.at(FElocPres)->getMapUnique());
1167 A->getBlock(0,1)->fillComplete(domainVec_.at(FElocPres)->getMapUnique(),domainVec_.at(FElocVel)->getMapVecFieldUnique());
1168 }
1169
1170 if(assembleMode == "Rhs"){
1171
1172 MultiVectorPtr_Type resVecUnique_u = Teuchos::rcp( new MultiVector_Type( domainVec_.at(FElocVel)->getMapVecFieldUnique(), 1 ) );
1173 MultiVectorPtr_Type resVecUnique_p = Teuchos::rcp( new MultiVector_Type( domainVec_.at(FElocPres)->getMapUnique(), 1 ) );
1174
1175 resVecUnique_u->putScalar(0.);
1176 resVecUnique_p->putScalar(0.);
1177
1178 resVecUnique_u->exportFromVector( resVec_u, true, "Add" );
1179 resVecUnique_p->exportFromVector( resVec_p, true, "Add" );
1180
1181 resVec->addBlock(resVecUnique_u,0);
1182 resVec->addBlock(resVecUnique_p,1);
1183 }
1184
1185
1186}
1187
1188
1193template <class SC, class LO, class GO, class NO>
1194void FE_ElementAssembly<SC,LO,GO,NO>::changeLinearizationFE(std::string linearization)
1195{
1196 for (UN T=0; T<assemblyFEElements_.size(); T++) // For each assembledFEElement change Linearization
1197 {
1198 assemblyFEElements_[T]->changeLinearization(linearization);
1199 }
1200}
1201
1202
1211
1212
1213template <class SC, class LO, class GO, class NO>
1214void FE_ElementAssembly<SC,LO,GO,NO>::computeSteadyViscosityFE_CM(int dim,
1215 std::string FETypeVelocity,
1216 std::string FETypePressure,
1217 int dofsVelocity,
1218 int dofsPressure,
1219 MultiVectorPtr_Type u_rep,
1220 MultiVectorPtr_Type p_rep,
1221 ParameterListPtr_Type params){
1222
1223
1224 UN FElocVel = checkFE(dim,FETypeVelocity); // Checks for different domains which belongs to a certain fetype
1225 UN FElocPres = checkFE(dim,FETypePressure); // Checks for different domains which belongs to a certain fetype
1226
1227 ElementsPtr_Type elements = domainVec_.at(FElocVel)->getElementsC();
1228 ElementsPtr_Type elementsPres = domainVec_.at(FElocPres)->getElementsC();
1229
1230 vec_dbl_Type solution(0);
1231 vec_dbl_Type solution_u;
1232 vec_dbl_Type solution_p;
1233 vec_dbl_Type solution_viscosity;
1234
1235 // We have to compute viscosity solution in each element
1236 MultiVectorPtr_Type Sol_viscosity = Teuchos::rcp( new MultiVector_Type( domainVec_.at(FElocVel)->getElementMap(), 1 ) ); //
1237 BlockMultiVectorPtr_Type visco_output = Teuchos::rcp( new BlockMultiVector_Type(1) );
1238 visco_output->addBlock(Sol_viscosity,0);
1239
1240
1241 for (UN T=0; T<assemblyFEElements_.size(); T++) {
1242
1243 vec_dbl_Type solution(0);
1244 solution_u = getSolution(elements->getElement(T).getVectorNodeList(), u_rep,dofsVelocity); // get the solution inside an element on the nodes
1245 solution_p = getSolution(elementsPres->getElement(T).getVectorNodeList(), p_rep,dofsPressure);
1246
1247 solution.insert( solution.end(), solution_u.begin(), solution_u.end() ); // here we insert the solution
1248 solution.insert( solution.end(), solution_p.begin(), solution_p.end() );
1249
1250 assemblyFEElements_[T]->updateSolution(solution); // here we update the value of the solutions inside an element
1251
1252 assemblyFEElements_[T]->computeLocalconstOutputField(); // we compute the viscosity inside an element
1253 solution_viscosity = assemblyFEElements_[T]->getLocalconstOutputField();
1254
1255 Teuchos::ArrayRCP<SC> resArray_block = visco_output->getBlockNonConst(0)->getDataNonConst(0); // First
1256 resArray_block[T] = solution_viscosity[0]; // although it is a vector it only has one entry because we compute the value in the center of the element
1257
1258 } // end loop over all elements
1259 // We could also instead of just overwrite it add an aditional block such that we could also compute other output fields and save it in there
1260 this->const_output_fields= visco_output;
1261
1262
1263}
1264
1265
1266
1277
1278template <class SC, class LO, class GO, class NO>
1279void FE_ElementAssembly<SC,LO,GO,NO>::addFeBlockMv(BlockMultiVectorPtr_Type &res, vec_dbl_ptr_Type rhsVec, FiniteElement elementBlock1,FiniteElement elementBlock2, int dofs1, int dofs2 ){
1280
1281 Teuchos::ArrayRCP<SC> resArray_block1 = res->getBlockNonConst(0)->getDataNonConst(0);
1282
1283 Teuchos::ArrayRCP<SC> resArray_block2 = res->getBlockNonConst(1)->getDataNonConst(0);
1284
1285 vec_LO_Type nodeList_block1 = elementBlock1.getVectorNodeList();
1286
1287 vec_LO_Type nodeList_block2 = elementBlock2.getVectorNodeList();
1288
1289 for(int i=0; i< nodeList_block1.size() ; i++){
1290 for(int d=0; d<dofs1; d++){
1291 resArray_block1[nodeList_block1[i]*dofs1+d] += (*rhsVec)[i*dofs1+d];
1292 }
1293 }
1294 int offset = nodeList_block1.size()*dofs1;
1295
1296 for(int i=0; i < nodeList_block2.size(); i++){
1297 for(int d=0; d<dofs2; d++)
1298 resArray_block2[nodeList_block2[i]*dofs2+d] += (*rhsVec)[i*dofs2+d+offset];
1299 }
1300
1301}
1302
1303
1316template <class SC, class LO, class GO, class NO>
1317void FE_ElementAssembly<SC,LO,GO,NO>::addFeBlock(BlockMatrixPtr_Type &A, SmallMatrixPtr_Type elementMatrix, FiniteElement element, MapConstPtr_Type mapRow, int row, int column, tuple_disk_vec_ptr_Type problemDisk){
1318
1319 int dofs1 = std::get<2>(problemDisk->at(row));
1320
1321 int numNodes1 = std::get<3>(problemDisk->at(row));
1322
1323 int dofsBlock1 = dofs1*numNodes1;
1324
1325 Teuchos::Array<SC> value( numNodes1, 0. );
1326 Teuchos::Array<GO> columnIndices( numNodes1, 0 );
1327
1328 for (UN i=0; i < numNodes1 ; i++) {
1329 for(int di=0; di<dofs1; di++){
1330 GO rowID =GO (dofs1* mapRow->getGlobalElement( element.getNode(i) )+di);
1331 for(int d=0; d<dofs1; d++){
1332 for (UN j=0; j < columnIndices.size(); j++){
1333 columnIndices[j] = GO ( dofs1 * mapRow->getGlobalElement( element.getNode(j) ) + d );
1334 value[j] = (*elementMatrix)[dofs1*i+di][dofs1*j+d];
1335 }
1336 A->getBlock(row,column)->insertGlobalValues( rowID, columnIndices(), value() ); // Automatically adds entries if a value already exists
1337 }
1338 }
1339 }
1340}
1341
1352template <class SC, class LO, class GO, class NO>
1353void FE_ElementAssembly<SC,LO,GO,NO>::initAssembleFEElements(std::string elementType,tuple_disk_vec_ptr_Type problemDisk,ElementsPtr_Type elements, ParameterListPtr_Type params,vec2D_dbl_ptr_Type pointsRep, MapConstPtr_Type elementMap){
1354
1355 vec2D_dbl_Type nodes;
1356 for (UN T=0; T<elements->numberElements(); T++) {
1357
1358 nodes = getCoordinates(elements->getElement(T).getVectorNodeList(), pointsRep);
1359
1360 AssembleFEFactory<SC,LO,GO,NO> assembleFEFactory;
1361
1362 AssembleFEPtr_Type assemblyFE = assembleFEFactory.build(elementType,elements->getElement(T).getFlag(),nodes, params,problemDisk);
1363 //
1364 assemblyFE->setGlobalElementID(elementMap->getGlobalElement(T));
1365
1366 assemblyFEElements_.push_back(assemblyFE);
1367
1368 }
1369
1370}
1371
1381
1382template <class SC, class LO, class GO, class NO>
1383vec2D_dbl_Type FE_ElementAssembly<SC,LO,GO,NO>::getCoordinates(vec_LO_Type localIDs, vec2D_dbl_ptr_Type points){
1384
1385 vec2D_dbl_Type coordinates(0,vec_dbl_Type( points->at(0).size()));
1386 for(int i=0; i < localIDs.size() ; i++){
1387 coordinates.push_back(points->at(localIDs[i]));
1388 }
1389
1390 return coordinates;
1391}
1392
1402
1403template <class SC, class LO, class GO, class NO>
1404vec_dbl_Type FE_ElementAssembly<SC,LO,GO,NO>::getSolution(vec_LO_Type localIDs, MultiVectorPtr_Type u_rep, int dofsVelocity){
1405
1406 Teuchos::ArrayRCP<SC> uArray = u_rep->getDataNonConst(0);
1407
1408 vec_dbl_Type solution(0);
1409 for(int i=0; i < localIDs.size() ; i++){
1410 for(int d=0; d<dofsVelocity; d++)
1411 solution.push_back(uArray[localIDs[i]*dofsVelocity+d]);
1412 }
1413
1414 return solution;
1415}
1416
1417
1418
1419template <class SC, class LO, class GO, class NO>
1420void FE_ElementAssembly<SC,LO,GO,NO>::assemblyEmptyMatrix(MatrixPtr_Type &A){
1421 A->fillComplete();
1422}
1423
1424
1425
1435template <class SC, class LO, class GO, class NO>
1436void FE_ElementAssembly<SC,LO,GO,NO>::assemblyLaplaceAssFE(int dim,
1437 std::string FEType,
1438 int degree,
1439 int dofs,
1440 BlockMatrixPtr_Type &A,
1441 bool callFillComplete,
1442 int FELocExternal){
1443 ParameterListPtr_Type params = Teuchos::getParametersFromXmlFile("parametersProblemLaplace.xml");
1444
1445 UN FEloc = checkFE(dim,FEType);
1446 ElementsPtr_Type elements = domainVec_.at(FEloc)->getElementsC();
1447 vec2D_dbl_ptr_Type pointsRep = domainVec_.at(FEloc)->getPointsRepeated();
1448 MapConstPtr_Type map = domainVec_.at(FEloc)->getMapRepeated();
1449 vec2D_dbl_Type nodes;
1450 int numNodes=dim+1;
1451 if(FEType == "P2"){
1452 numNodes= 6;
1453 if(dim==3)
1454 numNodes=10;
1455 }
1456 tuple_disk_vec_ptr_Type problemDisk = Teuchos::rcp(new tuple_disk_vec_Type(0));
1457 tuple_ssii_Type vel ("Laplace",FEType,dofs,numNodes);
1458 problemDisk->push_back(vel);
1459 if(assemblyFEElements_.size()== 0)
1460 initAssembleFEElements("Laplace",problemDisk,elements, params,pointsRep,domainVec_.at(0)->getElementMap());
1461 else if(assemblyFEElements_.size() != elements->numberElements())
1462 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Number Elements not the same as number assembleFE elements." );
1463 for (UN T=0; T<elements->numberElements(); T++) {
1464 assemblyFEElements_[T]->assembleJacobian();
1465 SmallMatrixPtr_Type elementMatrix = assemblyFEElements_[T]->getJacobian();
1466 addFeBlock(A, elementMatrix, elements->getElement(T), map, 0, 0, problemDisk);
1467
1468 }
1469 if(callFillComplete)
1470 A->getBlock(0,0)->fillComplete();
1471}
1472
1473
1474
1490
1491template <class SC, class LO, class GO, class NO>
1492void FE_ElementAssembly<SC,LO,GO,NO>::assemblyNonlinearLaplace(
1493 int dim, std::string FEType, int degree, MultiVectorPtr_Type u_rep,
1494 BlockMatrixPtr_Type &A, BlockMultiVectorPtr_Type &resVec,
1495 ParameterListPtr_Type params, std::string assembleMode, bool callFillComplete,
1496 int FELocExternal) {
1497
1498 ElementsPtr_Type elements = this->domainVec_.at(0)->getElementsC();
1499
1500 // Only scalar laplace
1501 int dofs = 1;
1502
1503 vec2D_dbl_ptr_Type pointsRep = this->domainVec_.at(0)->getPointsRepeated();
1504 MapConstPtr_Type map = this->domainVec_.at(0)->getMapRepeated();
1505
1506 vec_dbl_Type solution_u;
1507 vec_dbl_ptr_Type rhsVec;
1508
1509 int numNodes = 3;
1510 if (FEType == "P2") {
1511 numNodes = 6;
1512 }
1513 if (dim == 3) {
1514 numNodes = 4;
1515 if (FEType == "P2") {
1516 numNodes = 10;
1517 }
1518 }
1519
1520 // Tupel construction follows follwing pattern:
1521 // std::string: Physical Entity (i.e. Velocity) , std::string: Discretisation (i.e.
1522 // "P2"), int: Degrees of Freedom per Node, int: Number of Nodes per
1523 // element)
1524 tuple_disk_vec_ptr_Type problemDisk =
1525 Teuchos::rcp(new tuple_disk_vec_Type(0));
1526 tuple_ssii_Type temp("Solution", FEType, dofs, numNodes);
1527 problemDisk->push_back(temp);
1528
1529 // Construct an assembler for each element if not already done
1530 if (assemblyFEElements_.size() == 0) {
1531 initAssembleFEElements("NonLinearLaplace", problemDisk, elements, params, pointsRep, domainVec_.at(0)->getElementMap());
1532 } else if (assemblyFEElements_.size() != elements->numberElements()) {
1533 TEUCHOS_TEST_FOR_EXCEPTION(
1534 true, std::logic_error,
1535 "Number Elements not the same as number assembleFE elements.");
1536 }
1537
1538 MultiVectorPtr_Type resVec_u;
1539 BlockMultiVectorPtr_Type resVecRep;
1540
1541 if (assembleMode != "Rhs") {
1542 // add new or overwrite existing block (0,0) of system matrix
1543 // This is done in specific problem class for most other problems
1544 // Placing it here instead as more fitting
1545 auto A_block_zero_zero = Teuchos::rcp(
1546 new Matrix_Type(this->domainVec_.at(0)->getMapUnique(), this->domainVec_.at(0)->getApproxEntriesPerRow()));
1547
1548 A->addBlock(A_block_zero_zero, 0, 0);
1549 } else {
1550 // Or same for the residual vector
1551 resVec_u = Teuchos::rcp(new MultiVector_Type(map, 1));
1552 resVecRep = Teuchos::rcp(new BlockMultiVector_Type(1));
1553 resVecRep->addBlock(resVec_u, 0);
1554 }
1555 // Call assembly routines on each element
1556 for (UN T = 0; T < assemblyFEElements_.size(); T++) {
1557 vec_dbl_Type solution(0);
1558
1559 // Update solution on the element
1560 solution_u = getSolution(elements->getElement(T).getVectorNodeList(),
1561 u_rep, dofs);
1562 solution.insert(solution.end(), solution_u.begin(), solution_u.end());
1563 assemblyFEElements_[T]->updateSolution(solution);
1564
1565 if (assembleMode == "Jacobian") {
1566 SmallMatrixPtr_Type elementMatrix;
1567 assemblyFEElements_[T]->assembleJacobian();
1568 elementMatrix = assemblyFEElements_[T]->getJacobian();
1569 // Insert (additive) the local element Jacobian into the global
1570 // matrix
1571 assemblyFEElements_[T]
1572 ->advanceNewtonStep(); // n genereal non linear solver step
1573 addFeBlock(A, elementMatrix, elements->getElement(T), map, 0, 0,
1574 problemDisk);
1575 }
1576
1577 if (assembleMode == "Rhs") {
1578 assemblyFEElements_[T]->assembleRHS();
1579 rhsVec = assemblyFEElements_[T]->getRHS();
1580 // Name RHS comes from solving linear systems
1581 // For nonlinear systems RHS synonymous to residual
1582 // Insert (additive) the updated residual into the global residual
1583 // vector
1584 addFeBlockMv(resVecRep, rhsVec, elements->getElement(T), dofs);
1585 }
1586 }
1587 if (callFillComplete && assembleMode != "Rhs") {
1588 // Signal that editing A has finished. This causes the entries of A to
1589 // be redistributed across the MPI ranks
1590 A->getBlock(0, 0)->fillComplete(domainVec_.at(0)->getMapUnique(),
1591 domainVec_.at(0)->getMapUnique());
1592 }
1593 if (assembleMode == "Rhs") {
1594 // Export from overlapping residual to unique residual
1595 MultiVectorPtr_Type resVecUnique = Teuchos::rcp(
1596 new MultiVector_Type(domainVec_.at(0)->getMapUnique(), 1));
1597 resVecUnique->putScalar(0.);
1598 resVecUnique->exportFromVector(resVec_u, true, "Add");
1599 resVec->addBlock(resVecUnique, 0);
1600 }
1601}
1602
1603
1607template <class SC, class LO, class GO, class NO>
1608int FE_ElementAssembly<SC,LO,GO,NO>::checkFE(int dim,
1609 std::string FEType){
1610
1611 int FEloc;
1612 std::vector<int> matches;
1613 for (int i = 0; i < domainVec_.size(); i++) {
1614 if (domainVec_.at(i)->getDimension() == dim)
1615 matches.push_back(i);
1616 }
1617
1618 bool found = false;
1619 for (int i = 0; i < matches.size();i++) {
1620 if (domainVec_.at( matches.at(i) )->getFEType() == FEType) {
1621 FEloc = matches.at(i);
1622 found = true;
1623 }
1624 }
1625
1626 TEUCHOS_TEST_FOR_EXCEPTION(!found, std::logic_error ,"Combination of dimenson(2/3) and FE Type(P1/P2) not defined yet. Use addFE(domain)");
1627
1628 return FEloc;
1629}
1630
1631}
1632#endif
FEDD::FE_ElementAssembly::getCoordinates
vec2D_dbl_Type getCoordinates(vec_LO_Type localIDs, vec2D_dbl_ptr_Type points)
Returns coordinates of local node ids.
Definition FE_ElementAssembly_def.hpp:1383
FEDD::FE_ElementAssembly::checkFE
int checkFE(int Dimension, std::string FEType)
Checks which domain corresponds to certain FE Type and dimension.
Definition FE_ElementAssembly_def.hpp:1608
FEDD::FE_ElementAssembly::computeSteadyViscosityFE_CM
void computeSteadyViscosityFE_CM(int dim, std::string FETypeVelocity, std::string FETypePressure, int dofsVelocity, int dofsPressure, MultiVectorPtr_Type u_rep, MultiVectorPtr_Type p_rep, ParameterListPtr_Type params)
Postprocessing: Using a converged velocity solution -> compute averaged viscosity inside an element a...
Definition FE_ElementAssembly_def.hpp:1214
FEDD::FE_ElementAssembly::changeLinearizationFE
void changeLinearizationFE(std::string linearization)
Method to loop over all assembleFESpecific elements and set the defined linearization.
Definition FE_ElementAssembly_def.hpp:1194
FEDD::FE_ElementAssembly::getSolution
vec_dbl_Type getSolution(vec_LO_Type localIDs, MultiVectorPtr_Type u_rep, int dofsVelocity)
Returns entries of u of element.
Definition FE_ElementAssembly_def.hpp:1404
FEDD::FE_ElementAssembly::assemblyNonLinearElasticity
void assemblyNonLinearElasticity(int dim, std::string FEType, int degree, int dofs, MultiVectorPtr_Type d_rep, BlockMatrixPtr_Type &A, BlockMultiVectorPtr_Type &resVec, ParameterListPtr_Type params, bool callFillComplete=true, int FELocExternal=-1)
Assembly of Jacobian for nonlinear Elasticity.
Definition FE_ElementAssembly_def.hpp:398
FEDD::FE_ElementAssembly::assemblyNonlinearLaplace
void assemblyNonlinearLaplace(int dim, std::string FEType, int degree, MultiVectorPtr_Type u_rep, BlockMatrixPtr_Type &A, BlockMultiVectorPtr_Type &resVec, ParameterListPtr_Type params, std::string assembleMode, bool callFillComplete=true, int FELocExternal=-1)
Assembly of Jacobian for nonlinear Laplace example.
Definition FE_ElementAssembly_def.hpp:1492
FEDD::FE_ElementAssembly::assemblyLinearElasticity
void assemblyLinearElasticity(int dim, std::string FEType, int degree, int dofs, MultiVectorPtr_Type d_rep, BlockMatrixPtr_Type &A, BlockMultiVectorPtr_Type &resVec, ParameterListPtr_Type params, bool reAssemble, std::string assembleMode, bool callFillComplete=true, int FELocExternal=-1)
Assembly of Jacobian.
Definition FE_ElementAssembly_def.hpp:311
FEDD::FE_ElementAssembly::assemblyAceDeformDiffuBlock
void assemblyAceDeformDiffuBlock(int dim, std::string FETypeChem, std::string FETypeSolid, int degree, int dofsChem, int dofsSolid, MultiVectorPtr_Type c_rep, MultiVectorPtr_Type d_rep, BlockMatrixPtr_Type &A, int blockRow, int blockCol, BlockMultiVectorPtr_Type &resVec, int block, ParameterListPtr_Type params, std::string assembleMode, bool callFillComplete=true, int FELocExternal=-1)
Definition FE_ElementAssembly_def.hpp:770
FEDD::FE_ElementAssembly::assemblyLaplaceAssFE
void assemblyLaplaceAssFE(int dim, std::string FEType, int degree, int dofs, BlockMatrixPtr_Type &A, bool callFillComplete, int FELocExternal=-1)
Assembly of constant stiffness matix for laplacian operator .
Definition FE_ElementAssembly_def.hpp:1436
FEDD::FE_ElementAssembly::assemblyNavierStokes
void assemblyNavierStokes(int dim, std::string FETypeVelocity, std::string FETypePressure, int degree, int dofsVelocity, int dofsPressure, MultiVectorPtr_Type u_rep, MultiVectorPtr_Type p_rep, BlockMatrixPtr_Type &A, BlockMultiVectorPtr_Type &resVec, SmallMatrix_Type coeff, ParameterListPtr_Type params, bool reAssemble, std::string assembleMode, bool callFillComplete=true, int FELocExternal=-1)
Assembly of Jacobian for NavierStokes.
Definition FE_ElementAssembly_def.hpp:1026
FEDD::FE_ElementAssembly::assemblyAceDeformDiffu
void assemblyAceDeformDiffu(int dim, std::string FETypeChem, std::string FETypeSolid, int degree, int dofsChem, int dofsSolid, MultiVectorPtr_Type c_rep, MultiVectorPtr_Type d_rep, BlockMatrixPtr_Type &A, BlockMultiVectorPtr_Type &resVec, ParameterListPtr_Type params, std::string assembleMode, bool callFillComplete=true, int FELocExternal=-1)
Definition FE_ElementAssembly_def.hpp:592
FEDD::FiniteElement
Definition FiniteElement.hpp:17
FEDD::SmallMatrix
This class represents a templated small Matrix of type T. Primarily created for 2x2 and 3x3 matrices....
Definition SmallMatrix.hpp:20
FEDD
Adaptive Mesh Refinement.
Definition AdaptiveMeshRefinement_decl.hpp:36