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