Finite Element Domain Decomposition Library
FEDDLib
Loading...
Searching...
No Matches
AssembleFE_SCI_SMC_MLCK_def.hpp
1#ifndef AssembleFE_SCI_SMC_MLCK_DEF_hpp
2#define AssembleFE_SCI_SMC_MLCK_DEF_hpp
3
4#include "AssembleFE_SCI_SMC_MLCK_decl.hpp"
5
6#ifdef FEDD_HAVE_ACEGENINTERFACE
7#include <aceinterface.hpp>
8#endif
9
10#include <vector>
11//#include <iostream>
12
13namespace FEDD {
14
15template <class SC, class LO, class GO, class NO>
16AssembleFE_SCI_SMC_MLCK<SC,LO,GO,NO>::AssembleFE_SCI_SMC_MLCK(int flag, vec2D_dbl_Type nodesRefConfig, ParameterListPtr_Type params,tuple_disk_vec_ptr_Type tuple):
17AssembleFE<SC,LO,GO,NO>(flag, nodesRefConfig, params, tuple)
18{
19 /*
20 fA -Fibre angle_1 30e0,
21 $[Lambda]$C50 -LambdaC50_2 0.12e1
22 $[Gamma]$3 -Gamma3_3 0.9e0,
23 $[Lambda]$BarCDotMax -LambdaBarCDotMax_4 0.3387e-1
24 $[Lambda]$BarCDotMin -LambdaBarCDotMin_5 -0.3387e-1
25 $[Gamma]$2 -Gamma2_6 50e0,
26 $[Gamma]$1 -Gamma1_7 0.50247e0,
27 $[Eta]$1 -Eta1_8 0.18745e0,
28 Ca50 -Ca50_9 0.4e0
29 k2 -K2_10 0.2e0
30 k5 -K5_11 0.2e0,
31 k3 -K3_12 0.134e0
32 k4 -K4_13 0.166e-2
33 k7 -K7_14 0.66e-4
34 $[Kappa]$C -KappaC_15 0.14636600000000002e3
35 $[Beta]$1 -Beta1_16 0.10097e-2
36 $[Mu]$a -MuA_17 0.9291e1
37 $[Alpha]$ -Alpha_18 0.2668e2
38 $[Epsilon]$1 -Epsilon1_19 0.15173775e3
39 $[Epsilon]$2 -Epsilon2_20 0.27566199999999996e1
40 c1 -C1_21 0.1152507e2
41 $[Alpha]$1 -Alpha1_22 0.127631e1
42 $[Alpha]$2 -Alpha2_23 0.308798e1
43 p1 -P1_24 0.3e0
44 p3 -P3_25 0.2e0
45 c50 -C50_26 0.5e0
46 d0 -D0_27 0.6e-4
47 m -M_28 0e0
48 startTime -StartTime_29 1000e0
49 $[Rho]$0 -Density_30 1e0
50
51 */
52
53 fA_= this->params_->sublist("Parameter Solid").get("FA",30.e0); // ??
54 lambdaC50_ = this->params_->sublist("Parameter Solid").get("LambdaC50",0.12e1); // ??
55 gamma3_= this->params_->sublist("Parameter Solid").get("Gamma3",0.9e0);
56 lambdaBarCDotMax_= this->params_->sublist("Parameter Solid").get("LambdaBarCDotMax",0.3387e-1); // ??
57 lambdaBarCDotMin_= this->params_->sublist("Parameter Solid").get("LambdaBarCDotMin",-0.3387e-1); // ??
58 gamma2_ = this->params_->sublist("Parameter Solid").get("Gamma2",50.0e0); // ??
59 gamma1_ = this->params_->sublist("Parameter Solid").get("Gamma1",0.50247e0);
60 eta1_ = this->params_->sublist("Parameter Solid").get("Eta1",0.18745e0); // ??
61 ca50_ = this->params_->sublist("Parameter Solid").get("Ca50",0.4e0); // ??
62 k2_ = this->params_->sublist("Parameter Solid").get("K2",0.2e0);
63 k5_ = this->params_->sublist("Parameter Solid").get("K5",0.2e0);
64 k3_ = this->params_->sublist("Parameter Solid").get("K3",0.134e0); // ??
65 k4_ = this->params_->sublist("Parameter Solid").get("K4",0.166e-2); // ??
66 k7_= this->params_->sublist("Parameter Solid").get("K7",0.66e-4); // ??
67 kappaC_ = this->params_->sublist("Parameter Solid").get("KappaC",146.36600000000002e0);
68 beta1_ = this->params_->sublist("Parameter Solid").get("Beta1",0.10097e-2); // ??
69 muA_ = this->params_->sublist("Parameter Solid").get("MuA",0.9291e1);
70 alpha_ = this->params_->sublist("Parameter Solid").get("Alpha",0.2668e2);
71 epsilon1_ = this->params_->sublist("Parameter Solid").get("Epsilon1", 0.15173775e3);
72 epsilon2_ = this->params_->sublist("Parameter Solid").get("Epsilon2",0.27566199999999996e1); // ??
73 c1_ = this->params_->sublist("Parameter Solid").get("C1",11.52507e0);
74 alpha1_ = this->params_->sublist("Parameter Solid").get("Alpha1",1.27631e0);
75 alpha2_ = this->params_->sublist("Parameter Solid").get("Alpha2",0.308798e1); // ??
76 p1_ = this->params_->sublist("Parameter Solid").get("P1",0.3e0);
77 p3_ = this->params_->sublist("Parameter Solid").get("P3",0.2e0);
78 c50_ = this->params_->sublist("Parameter Solid").get("C50",0.5e0);
79 d0_ = this->params_->sublist("Parameter Diffusion").get("D0",6.e-05);
80 m_ = this->params_->sublist("Parameter Solid").get("m",0.e0);
81 startTime_ = this->params_->sublist("Parameter Solid").get("ActiveStartTime",1001.e0); // At Starttime 1000 the diffused drug influences the material model. -> Active response at T=starttime
82 rho_ = this->params_->sublist("Parameter Solid").get("Rho",1.e0);
83
84 // iCode_ = this->params_->sublist("Parameter Solid").get("Intergration Code",18);
85 iCode_=18; //Only works for 18 currently!!
86
87 FEType_ = std::get<1>(this->diskTuple_->at(0)); // FEType of Disk
88 dofsSolid_ = std::get<2>(this->diskTuple_->at(0)); // Degrees of freedom per node
89 dofsChem_ = std::get<2>(this->diskTuple_->at(1)); // Degrees of freedom per node
90
91 numNodesSolid_ = std::get<3>(this->diskTuple_->at(0)); // Number of nodes of element
92 numNodesChem_ = std::get<3>(this->diskTuple_->at(1)); // Number of nodes of element
93
94 dofsElement_ = dofsSolid_*numNodesSolid_ + dofsChem_*numNodesChem_; // "Dimension of return matrix"
95
96 // Einlesen durch Parameterdatei irgendwann cool
97 history_ ={1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0}; // 48 values, 12 variables, 4 gausspoints
98 //.resize(48);
99
100 historyUpdated_.resize(48,0.);
101
102 solutionC_n_.resize(10,0.);
103 solutionC_n1_.resize(10,0.);
104
105 this->solution_.reset( new vec_dbl_Type ( dofsElement_,0.) );
106
107 /*timeParametersVec_.resize(0, vec_dbl_Type(2));
108 numSegments_ = this->params_->sublist("Timestepping Parameter").sublist("Timestepping Intervalls").get("Number of Segments",0);
109
110 for(int i=1; i <= numSegments_; i++){
111
112 double startTime = this->params_->sublist("Timestepping Parameter").sublist("Timestepping Intervalls").sublist(std::to_string(i)).get("Start Time",0.);
113 double dtTmp = this->params_->sublist("Timestepping Parameter").sublist("Timestepping Intervalls").sublist(std::to_string(i)).get("dt",0.1);
114
115 vec_dbl_Type segment = {startTime,dtTmp};
116 timeParametersVec_.push_back(segment);
117 }*/
118
119}
120
121template <class SC, class LO, class GO, class NO>
122void AssembleFE_SCI_SMC_MLCK<SC,LO,GO,NO>::assembleJacobian() {
123
124 SmallMatrixPtr_Type elementMatrix = Teuchos::rcp( new SmallMatrix_Type(dofsElement_,0.));
125
126 assemble_SCI_SMC_MLCK(elementMatrix);
127
128 this->jacobian_ = elementMatrix;
129
130}
131template <class SC, class LO, class GO, class NO>
132void AssembleFE_SCI_SMC_MLCK<SC,LO,GO,NO>::advanceInTime( double dt){
133
134 // If we have a time segment setting we switch to the demanded time increment
135 /*for(int i=0; i<numSegments_ ; i++){
136 if(this->timeStep_ +1.0e-12 > timeParametersVec_[i][0])
137 this->timeIncrement_=timeParametersVec_[i][1];
138 }*/
139 this->timeStep_ = this->timeStep_ + this->timeIncrement_;
140
141 this->timeIncrement_ = dt;
142
143 for(int i=0; i< 48; i++){
144 // if(this->timeStep_ > startTime_ +dt )
145 history_[i] = historyUpdated_[i];
146
147 }
148
149 for(int i=0; i< 10 ; i++)
150 solutionC_n_[i]=(*this->solution_)[i+30]; // this is the LAST solution of newton iterations
151
152
153}
154
155template <class SC, class LO, class GO, class NO>
156void AssembleFE_SCI_SMC_MLCK<SC,LO,GO,NO>::assembleRHS(){
157
158 this->rhsVec_.reset( new vec_dbl_Type ( dofsElement_,0.) );
159
160#ifdef FEDD_HAVE_ACEGENINTERFACE
161
162
163 double deltaT=this->getTimeIncrement();
164
165 double positions[30];
166 int count = 0;
167 //cout << "Positions " << endl;
168 for(int i=0;i<10;i++){
169 for(int j=0;j<3;j++){
170 positions[count] = this->getNodesRefConfig()[i][j];
171 count++;
172 // cout << " i " << count-1 << " " << positions[count-1] << endl;
173 }
174 }
175 //cout << endl;
176
177 double displacements[30];
178 for(int i = 0; i < 30; i++)
179 {
180 displacements[i]=(*this->solution_)[i];
181 }
182
183 double history[48];// = {1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0}; // 48 values, 12 variables, 4 gausspoints
184 for(int i = 0; i < history_.size(); i++)
185 history[i] = history_[i];
186
187
188 double concentrations[10];
189 for(int i = 0; i < 10; i++)
190 {
191 concentrations[i]= (*this->solution_)[i+30];
192 solutionC_n1_[i] = (*this->solution_)[i+30]; // in each newtonstep solution for n+1 is updated.
193 }
194
195 double accelerations[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
196
197 double domainData[] = {fA_, lambdaC50_, gamma3_, lambdaBarCDotMax_, lambdaBarCDotMin_,
198 gamma2_ , gamma1_, eta1_, ca50_, k2_, k5_, k3_, k4_, k7_ , kappaC_ ,
199 beta1_ , muA_ , alpha_, epsilon1_,epsilon2_,c1_,alpha1_,alpha2_,
200 p1_,p3_,c50_,d0_,m_,startTime_,rho_};
201
202 double rates[10];
203
204 for(int i=0; i<10 ; i++){
205 rates[i] =(solutionC_n1_[i]-solutionC_n_[i])/deltaT;//(solutionC_n1_[i])/deltaT; //-solutionC_n_[i](solutionC_n1_[i]-solutionC_n_[i])/deltaT;//
206 }
207
208
209 double time = this->getTimeStep()+deltaT;
210 double subIterationTolerance = 1.e-7;
211
212 // immer speicher und wenn es konvergiert, dann zur history machen
213 double historyUpdated[] = {1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0}; // 48 values, 12 variables, 4 gausspoints
214
215 AceGenInterface::DeformationDiffusionSmoothMuscleMLCKOnlyTetrahedra3D10 elem(positions, displacements, concentrations, accelerations, rates, domainData, history, subIterationTolerance, deltaT, time, iCode_);
216 elem.compute();
217
218 double *residuumRint = elem.getResiduumVectorRint();
219
220 double *residuumRDyn = elem.getResiduumVectorRdyn();
221
222 // getResiduumVectorRint(&positions[0], &displacements[0], &concentrations[0], &accelerations[0], &rates[0], &domainData[0], &history[0], subIterationTolerance, deltaT, time, iCode_, &historyUpdated[0], residuumRint);
223
224 // getResiduumVectorRdyn(&positions[0], &displacements[0], &concentrations[0], &accelerations[0],&rates[0], &domainData[0], &history[0], subIterationTolerance, deltaT, time, iCode_, &historyUpdated[0], residuumRDyn);
225 for(int i=0; i< 30 ; i++){
226 (*this->rhsVec_)[i] = residuumRint[i]; //+residuumRDyn[i];
227 }
228 double *residuumRc = elem.getResiduumVectorRc();
229 // getResiduumVectorRc(&positions[0], &displacements[0], &concentrations[0], &accelerations[0], &rates[0], &domainData[0], &history[0], subIterationTolerance, deltaT, time, iCode_, &historyUpdated[0], residuumRc);
230
231
232 for(int i=0; i< 10 ; i++){
233 (*this->rhsVec_)[i+30] = residuumRc[i];
234 }
235
236 // free(residuumRc);
237 // free(residuumRint);
238 // free(residuumRDyn);
239
240
241
242#endif
243
244}
245
262
263
264template <class SC,class LO, class GO, class NO>
265void AssembleFE_SCI_SMC_MLCK<SC,LO,GO,NO>::assemble_SCI_SMC_MLCK(SmallMatrixPtr_Type &elementMatrix){
266 // double deltat=this->getTimeIncrement();
267 // std::vector<double> deltat(1);
268 // deltat[0]=this->getTimeIncrement();
269
270#ifdef FEDD_HAVE_ACEGENINTERFACE
271
272 double deltaT=this->getTimeIncrement();
273
274 double positions[30];
275 int count = 0;
276 for(int i=0;i<10;i++){
277 for(int j=0;j<3;j++){
278 positions[count] = this->getNodesRefConfig()[i][j];
279 count++;
280 }
281 }
282 double displacements[30];
283 for(int i = 0; i < 30; i++)
284 {
285 displacements[i]= (*this->solution_)[i];
286
287 }
288 double history[48];// = {1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0}; // 48 values, 12 variables, 4 gausspoints
289 for(int i = 0; i < history_.size(); i++)
290 history[i] = history_[i]; //if (integrationCode==19)
291
292
293 double concentrations[10];
294 for(int i = 0; i < 10; i++)
295 {
296 concentrations[i]=(*this->solution_)[i+30];
297 solutionC_n1_[i]=(*this->solution_)[i+30]; // in each newtonstep solution for n+1 is updated.
298 }
299
300 double accelerations[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
301
302 double domainData[] = {fA_, lambdaC50_, gamma3_, lambdaBarCDotMax_, lambdaBarCDotMin_,
303 gamma2_ , gamma1_, eta1_, ca50_, k2_, k5_, k3_, k4_, k7_ , kappaC_ ,
304 beta1_ , muA_ , alpha_, epsilon1_,epsilon2_,c1_,alpha1_,alpha2_,
305 p1_,p3_,c50_,d0_,m_,startTime_,rho_};
306
307 double rates[10];
308 for(int i=0; i<10 ; i++){
309 rates[i] =(solutionC_n1_[i]-solutionC_n_[i]) / deltaT;//
310 }
311 // ##########################
312
313 double time = this->getTimeStep()+deltaT;
314 double subIterationTolerance = 1.e-7;
315
316 // immer speicher und wenn es konvergiert, dann zur history machen
317 // double historyUpdated[] = {1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0}; // 48 values, 12 variables, 4 gausspoints
318
319 AceGenInterface::DeformationDiffusionSmoothMuscleMLCKOnlyTetrahedra3D10 elem(positions, displacements, concentrations, accelerations, rates, domainData, history, subIterationTolerance, deltaT, time, iCode_);
320 elem.compute();
321
322 double *historyUpdated = elem.getHistoryUpdated();
323
324 double** stiffnessMatrixKuu = elem.getStiffnessMatrixKuu();
325 double** stiffnessMatrixKuc = elem.getStiffnessMatrixKuc();
326 double** stiffnessMatrixKcu = elem.getStiffnessMatrixKcu();
327 double** stiffnessMatrixKcc = elem.getStiffnessMatrixKcc();
328 double** massMatrixMc = elem.getMassMatrixMc();
329
330
331 for(int i=0; i< 48; i++){
332 historyUpdated_[i] = historyUpdated[i];
333 }
334
335 for(int i=0; i< 30; i++){
336 for(int j=0; j<30; j++){
337 //if(std::fabs(stiffnessMatrixKuu[i][j]) > 1e7)
338 // cout << " !!! Sus entry Kuu [" << i << "][" << j << "] " << stiffnessMatrixKuu[i][j] << endl;
339
340 (*elementMatrix)[i][j]=stiffnessMatrixKuu[i][j];
341 }
342 }
343 for(int i=0; i< 30; i++){
344 for(int j=0; j<10; j++){
345 //if(std::fabs(stiffnessMatrixKuc[i][j]) > 1e7)
346 // cout << " !!! Sus entry Kuc [" << i << "][" << j << "] " << stiffnessMatrixKuc[i][j] << endl;
347
348 (*elementMatrix)[i][j+30]=stiffnessMatrixKuc[i][j];
349 }
350 }
351 for(int i=0; i< 10; i++){
352 for(int j=0; j<30; j++){
353 //if(std::fabs(stiffnessMatrixKcu[i][j]) > 1e7)
354 // cout << " !!! Sus entry Kcu [" << i << "][" << j << "] " << stiffnessMatrixKcu[i][j] << endl;
355
356 (*elementMatrix)[i+30][j]=stiffnessMatrixKcu[i][j];
357 }
358 }
359 for(int i=0; i< 10; i++){
360 for(int j=0; j<10; j++){
361 //if(std::fabs(massMatrixMc[i][j]) > 1e5 || std::fabs(stiffnessMatrixKcc[i][j]) > 1e5 )
362 // cout << " !!! Sus entry Mass [" << i << "][" << j << "] " << massMatrixMc[i][j] << " or stiff Kcc " << stiffnessMatrixKcc[i][j] << endl;
363
364 (*elementMatrix)[i+30][j+30] =stiffnessMatrixKcc[i][j] +(1./deltaT)*massMatrixMc[i][j]; //
365 }
366 }
367
368
369
370
371#endif
372
373
374}
375
376
377} // namespace FEDD
378#endif // AssembleFE_SCI_SMC_MLCK_DEF_hpp
FEDD::AssembleFE_SCI_SMC_MLCK::assembleJacobian
void assembleJacobian() override
Assemble the element Jacobian matrix.
Definition AssembleFE_SCI_SMC_MLCK_def.hpp:122
FEDD::AssembleFE_SCI_SMC_MLCK::advanceInTime
void advanceInTime(double dt) override
This function is called every time the FEDDLib proceeds from one to the next time step....
Definition AssembleFE_SCI_SMC_MLCK_def.hpp:132
FEDD::AssembleFE_SCI_SMC_MLCK::assembleRHS
void assembleRHS() override
Assemble the element right hand side vector.
Definition AssembleFE_SCI_SMC_MLCK_def.hpp:156
FEDD::AssembleFE
This abstract class defining the interface for any type of element assembly rountines in the FEDDLib.
Definition AssembleFE_decl.hpp:61
FEDD::AssembleFE< default_sc, default_lo, default_go, default_no >::getTimeIncrement
double getTimeIncrement()
Definition AssembleFE_decl.hpp:204
FEDD::AssembleFE< default_sc, default_lo, default_go, default_no >::getNodesRefConfig
vec2D_dbl_Type getNodesRefConfig()
Definition AssembleFE_def.hpp:144
FEDD
Adaptive Mesh Refinement.
Definition AdaptiveMeshRefinement.cpp:5