Finite Element Domain Decomposition Library
FEDDLib
Loading...
Searching...
No Matches
AssembleFE_SCI_SMC_Active_Growth_Reorientation_def.hpp
1#ifndef AssembleFE_SCI_SMC_Active_Growth_Reorientation_DEF_hpp
2#define AssembleFE_SCI_SMC_Active_Growth_Reorientation_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_Active_Growth_Reorientation<SC,LO,GO,NO>::AssembleFE_SCI_SMC_Active_Growth_Reorientation(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 $[Gamma]$6 -Gamma6_24 0.15e1
42 $[Lambda]$P50 -LambdaP50_25 1.e0
43 kDotMin -KDotMin_26 -0.118863e-2
44 $[Zeta]$1 -Zeta1_27 100.e0
45 kDotMax -KDotMax_28 0.51028e-3
46 $[Gamma]$4 -Gamma4_29 200.e0
47 $[Lambda]$BarDotPMin -LambdaBarDotPMin_30 -0.17689e-3
48 $[Lambda]$BarDotPMax -LambdaBarDotPMax_31 0.13957e-3
49 $[Gamma]$5 -Gamma5_32 50.e0
50 $[Zeta]$2 -Zeta2_33 1000.e0
51 $[CapitalDelta]$$[Lambda]$BarPMin -DeltaLambdaBarPMin_34 -0.1e-4
52 p1 -P1_35 0.3e0
53 p3 -P3_36 0.2e0
54 c50 -C50_37 0.5e0
55 d0 -D0_38 0.6e-4
56 m -M_39 0e0
57 activeStartTime -ActiveStartTime_40 1000.e0
58 k$[Eta]$Plus -kEtaPlus_41 0.6e0
59 m$[Eta]$Plus -mEtaPlus_42 5.e0
60 growthStartTime -growthStartTime_43 1.e0
61 reorientationStartTime -reorientationStartTime_44 1.e0
62 growthEndTime -growthEndTime_45 100.e0
63 reorientationEndTime -reorientationEndTime_46 100.e0
64 k$[Theta]$Plus -KThetaPlus_47 1.e0
65 k$[Theta]$Minus -KThetaMinus_48 1.e0
66 m$[Theta]$Plus -MThetaPlus_49 3.e0
67 m$[Theta]$Minus -MThetaMinus_50 3.e0
68 $[Theta]$Plus1 -ThetaPlus1_51 0.1000882e1
69 $[Theta]$Plus2 -ThetaPlus2_52 0.1234826e1
70 $[Theta]$Plus3 -ThetaPlus3_53 0.11414189999999999e1
71 $[Theta]$Minus1 -ThetaMinus1_54 0.98e0
72 $[Theta]$Minus2 -ThetaMinus2_55 0.98e0
73 $[Theta]$Minus3 -ThetaMinus3_56 0.98e0
74 $[Rho]$ -Density_57 1.e0
75
76 */
77
78
79 // -------------------- Parameter ---------------------
80 fA_= this->params_->sublist("Parameter Solid").get("FA",30.e0); // ??
81 lambdaC50_ = this->params_->sublist("Parameter Solid").get("LambdaC50",0.12e1); // ??
82 gamma3_= this->params_->sublist("Parameter Solid").get("Gamma3",0.9e0);
83 lambdaBarCDotMax_= this->params_->sublist("Parameter Solid").get("LambdaBarCDotMax",0.443e-1); // ??
84 lambdaBarCDotMin_= this->params_->sublist("Parameter Solid").get("LambdaBarCDotMin",-0.443e-1); // ??
85 gamma2_ = this->params_->sublist("Parameter Solid").get("Gamma2",50.0e0); // ??
86 gamma1_ = this->params_->sublist("Parameter Solid").get("Gamma1",0.5131e0);
87 eta_ = this->params_->sublist("Parameter Solid").get("Eta",0.18745e0); // ??
88 ca50_ = this->params_->sublist("Parameter Solid").get("Ca50",0.4e0); // ??
89 k2_ = this->params_->sublist("Parameter Solid").get("K2",0.2e0);
90 k5_ = this->params_->sublist("Parameter Solid").get("K5",0.2e0);
91 k3_ = this->params_->sublist("Parameter Solid").get("K3",0.134e0); // ??
92 k4_ = this->params_->sublist("Parameter Solid").get("K4",0.166e-2); // ??
93 k7_= this->params_->sublist("Parameter Solid").get("K7",0.66e-4); // ??
94 kappa_ = this->params_->sublist("Parameter Solid").get("Kappa",0.148262e0);
95 beta1_ = this->params_->sublist("Parameter Solid").get("Beta1",0.1006e-2); // ??
96 muA_ = this->params_->sublist("Parameter Solid").get("MuA",0.11857e-1);
97 beta2_ = this->params_->sublist("Parameter Solid").get("Beta2",0.2668e-1);
98 alpha2_ = this->params_->sublist("Parameter Solid").get("Alpha2",0.15173775e0);
99 alpha3_ = this->params_->sublist("Parameter Solid").get("Alpha3",0.275662e1); // ??
100 alpha1_ = this->params_->sublist("Parameter Solid").get("Alpha1",11.52507e-3);
101 alpha4_ = this->params_->sublist("Parameter Solid").get("Alpha4",1.27631e-3);
102 alpha5_ = this->params_->sublist("Parameter Solid").get("Alpha5",0.308798e1); // ??
103 gamma6_ = this->params_->sublist("Parameter Solid").get("Gamma6",0.15e1);
104 lambdaP50_ = this->params_->sublist("Parameter Solid").get("LambdaP50",1.0e0);
105 kDotMin_ = this->params_->sublist("Parameter Solid").get("KDotMin",-0.10694e-1);
106 zeta1_ = this->params_->sublist("Parameter Solid").get("Zeta1",100.e0);
107 kDotMax_ = this->params_->sublist("Parameter Solid").get("KDotMax",0.9735e-3);
108 gamma4_ = this->params_->sublist("Parameter Solid").get("Gamma4",200.e0);
109 lambdaBarDotPMin_ = this->params_->sublist("Parameter Solid").get("LambdaBarDotMin",-0.2323e-3);
110 lambdaBarDotPMax_ = this->params_->sublist("Parameter Solid").get("LambdaBarDotMax",0.699e-4);
111 gamma5_ = this->params_->sublist("Parameter Solid").get("Gamma5", 50.e0 );
112 zeta2_ = this->params_->sublist("Parameter Solid").get("Zeta2", 1000.e0 );
113 DeltaLambdaBarPMin_ =this->params_->sublist("Parameter Solid").get("DeltaLambdaBarPMin", -0.1e-4);
114 p1_ = this->params_->sublist("Parameter Solid").get("P1",0.6e0);
115 p3_ = this->params_->sublist("Parameter Solid").get("P3",0.6e0);
116 c50_ = this->params_->sublist("Parameter Solid").get("C50",0.5e0);
117 d0_ = this->params_->sublist("Parameter Diffusion").get("D0",6.e-05);
118 m_ = this->params_->sublist("Parameter Solid").get("m",0.e0);
119 activeStartTime_ = this->params_->sublist("Parameter Solid").get("ActiveStartTime",1.0); // At Starttime 1000 the diffused drug influences the material model. -> Active response at T=starttime
120 kEtaPlus_ = this->params_->sublist("Parameter Solid").get("KEtaPlus",0.1e-3);
121 mEtaPlus_ = this->params_->sublist("Parameter Solid").get("MEtaPlus",5.0e0);
122 growthStartTime_ = this->params_->sublist("Parameter Solid").get("GrowthStartTime",0.e0);
123 reorientationStartTime_ = this->params_->sublist("Parameter Solid").get("ReorientationStartTime",0.e0);
124 growthEndTime_ = this->params_->sublist("Parameter Solid").get("GrowthEndTime",0.e0);
125 reorientationEndTime_ = this->params_->sublist("Parameter Solid").get("ReorientationEndTime",0.e0);
126 kThetaPlus_ = this->params_->sublist("Parameter Solid").get("KThetaPlus",1.e-4);
127 kThetaMinus_ = this->params_->sublist("Parameter Solid").get("KThetaMinus",1.e-4);
128 mThetaPlus_ = this->params_->sublist("Parameter Solid").get("MThetaPlus",3.e0);
129 mThetaMinus_ = this->params_->sublist("Parameter Solid").get("MThetaMinus",3.e0);
130 thetaPlus1_ = this->params_->sublist("Parameter Solid").get("ThetaPlus1",1.005063);
131 thetaPlus2_ = this->params_->sublist("Parameter Solid").get("ThetaPlus2",1.020595);
132 thetaPlus3_ = this->params_->sublist("Parameter Solid").get("ThetaPlus3",1.119918);
133 thetaMinus1_ = this->params_->sublist("Parameter Solid").get("ThetaMinus1",0.98e0);
134 thetaMinus2_ = this->params_->sublist("Parameter Solid").get("ThetaMinus2",0.98e0);
135 thetaMinus3_ = this->params_->sublist("Parameter Solid").get("ThetaMinus3",0.98e0);
136 kMin_ = this->params_->sublist("Parameter Solid").get("KMin",0.15e0);
137 rho_ = this->params_->sublist("Parameter Solid").get("Rho",1.e0);
138
139// iCode_ = this->params_->sublist("Parameter Solid").get("Intergration Code",18);
140 iCode_=18; //Only works for 18 currently!!
141
142 FEType_ = std::get<1>(this->diskTuple_->at(0)); // FEType of Disk
143 dofsSolid_ = std::get<2>(this->diskTuple_->at(0)); // Degrees of freedom per node
144 dofsChem_ = std::get<2>(this->diskTuple_->at(1)); // Degrees of freedom per node
145
146 numNodesSolid_ = std::get<3>(this->diskTuple_->at(0)); // Number of nodes of element
147 numNodesChem_ = std::get<3>(this->diskTuple_->at(1)); // Number of nodes of element
148
149 dofsElement_ = dofsSolid_*numNodesSolid_ + dofsChem_*numNodesChem_; // "Dimension of return matrix"
150
151 // Einlesen durch Parameterdatei irgendwann cool
152 this->history_ ={1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 1., 1., 1.82758, 1.82758, 1., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
153 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 1., 1., 1.82758, 1.82758, 1., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
154 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 1., 1., 1.82758, 1.82758, 1., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
155 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 1., 1., 1.82758, 1.82758, 1., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
156 /*{1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 1., 1., 1.512656, 1.512656, 1., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
157 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 1., 1., 1.512656, 1.512656, 1., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
158 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 1., 1., 1.512656, 1.512656, 1., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
159 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 1., 1., 1.512656, 1.512656, 1., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};*/
160
161 this->historyLength_ = 0;
162#ifdef FEDD_HAVE_ACEGENINTERFACE
163 AceGenInterface::DeformationDiffusionSmoothMuscleActiveGrowthReorientationTetrahedra3D10 tempElem(iCode_);
164 this->historyLength_ = tempElem.getHistoryLength();
165
166 TEUCHOS_TEST_FOR_EXCEPTION(this->historyLength_ != this->history_.size(), std::logic_error, "History input length does not match history size of model! \n Hisory input length: " << this->history_.size() << "\n History size of model: " << this->historyLength_ << "\n");
167#endif
168
169 this->historyUpdated_.resize(this->historyLength_,0.);
170
171
172 solutionC_n_.resize(10,0.);
173 solutionC_n1_.resize(10,0.);
174
175 this->solution_.reset( new vec_dbl_Type ( dofsElement_,0.) );
176
177 /*timeParametersVec_.resize(0, vec_dbl_Type(2));
178 numSegments_ = this->params_->sublist("Timestepping Parameter").sublist("Timestepping Intervalls").get("Number of Segments",0);
179
180 for(int i=1; i <= numSegments_; i++){
181
182 double startTime = this->params_->sublist("Timestepping Parameter").sublist("Timestepping Intervalls").sublist(std::to_string(i)).get("Start Time",0.);
183 double dtTmp = this->params_->sublist("Timestepping Parameter").sublist("Timestepping Intervalls").sublist(std::to_string(i)).get("dt",0.1);
184
185 vec_dbl_Type segment = {startTime,dtTmp};
186 timeParametersVec_.push_back(segment);
187 }*/
188 massMatrix_ = Teuchos::rcp( new SmallMatrix_Type(10,0.));
189}
190
191template <class SC, class LO, class GO, class NO>
192void AssembleFE_SCI_SMC_Active_Growth_Reorientation<SC,LO,GO,NO>::assembleJacobian() {
193
194 SmallMatrixPtr_Type elementMatrix = Teuchos::rcp( new SmallMatrix_Type(this->dofsElement_,0.));
195
196 assemble_SCI_SMC_Active_Growth_Reorientation(elementMatrix);
197
198 this->jacobian_ = elementMatrix;
199
200}
201template <class SC, class LO, class GO, class NO>
202void AssembleFE_SCI_SMC_Active_Growth_Reorientation<SC,LO,GO,NO>::advanceInTime( double dt){
203
204 // If we have a time segment setting we switch to the demanded time increment
205 /*for(int i=0; i<numSegments_ ; i++){
206 if(this->timeStep_+1.0e-12 > timeParametersVec_[i][0])
207 this->timeIncrement_=timeParametersVec_[i][1];
208 }*/
209
210 this->timeStep_ = this->timeStep_ + this->timeIncrement_;
211
212 this->timeIncrement_ = dt;
213
214 for(int i=0; i< this->historyLength_; i++){
215 //if(this->timeStep_ > activeStartTime_ +dt )
216 this->history_[i] = this->historyUpdated_[i];
217 }
218
219 for(int i=0; i< 10 ; i++)
220 this->solutionC_n_[i]=(*this->solution_)[i+30]; // this is the LAST solution of newton iterations
221
222
223}
224
225template <class SC, class LO, class GO, class NO>
226void AssembleFE_SCI_SMC_Active_Growth_Reorientation<SC,LO,GO,NO>::assembleRHS(){
227
228 this->rhsVec_.reset( new vec_dbl_Type ( this->dofsElement_,0.) );
229
230#ifdef FEDD_HAVE_ACEGENINTERFACE
231
232
233 double deltaT=this->getTimeIncrement();
234
235 double positions[30];
236 int count = 0;
237 //cout << "Positions " << endl;
238 for(int i=0;i<10;i++){
239 for(int j=0;j<3;j++){
240 positions[count] = this->getNodesRefConfig()[i][j];
241 count++;
242 // cout << " | " << positions[count-1] ;
243 }
244 // cout << endl;
245
246 }
247 //cout << " --- " << endl;
248 // std::ofstream myfile;
249 // myfile.open ("outputs.txt",ios::app);
250 // myfile << "Positions " << endl;
251 // for(int i=0;i<30;i++)
252 // myfile << positions[i] << endl;
253 // myfile << endl;
254
255
256 double displacements[30];
257 for(int i = 0; i < 30; i++)
258 {
259 displacements[i]=(*this->solution_)[i];
260 }
261
262
263 std::vector<double> history(this->historyLength_, 0.0);
264 for(int i = 0; i < this->historyLength_; i++){
265 history[i] = this->history_[i];
266 }
267
268 double concentrations[10];
269 for(int i = 0; i < 10; i++)
270 {
271 concentrations[i]= (*this->solution_)[i+30];
272 this->solutionC_n1_[i] = (*this->solution_)[i+30]; // in each newtonstep solution for n+1 is updated.
273 }
274
275 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};
276
277 std::vector<double> domainData = {this->fA_, this->lambdaC50_, this->gamma3_, this->lambdaBarCDotMax_, this->lambdaBarCDotMin_, this->gamma2_, this->gamma1_, this->eta_, this->ca50_, this->k2_, this->k5_,
278 this->k3_, this->k4_, this->k7_, this->kappa_, this->beta1_, this->muA_, this->beta2_, this->alpha2_, this->alpha3_, this->alpha1_, this->alpha4_, this->alpha5_,
279 this->gamma6_, this->lambdaP50_, this->kDotMin_, this->zeta1_, this->kDotMax_, this->gamma4_, this->lambdaBarDotPMin_, this->lambdaBarDotPMax_, this->gamma5_, this->zeta2_,
280 this->DeltaLambdaBarPMin_, this->p1_, this->p3_, this->c50_, this->d0_, this->m_, this->activeStartTime_, this->kEtaPlus_, this->mEtaPlus_, this->growthStartTime_,
281 this->reorientationStartTime_, this->growthEndTime_, this->reorientationEndTime_, this->kThetaPlus_, this->kThetaMinus_, this->mThetaPlus_, this->mThetaMinus_, this->thetaPlus1_,
282 this->thetaPlus2_, this->thetaPlus3_, this->thetaMinus1_, this->thetaMinus2_, this->thetaMinus3_, this->kMin_, this->rho_};
283
284 double rates[10];
285
286 for(int i=0; i<10 ; i++){
287 rates[i] =(this->solutionC_n1_[i]-this->solutionC_n_[i])/deltaT;//(solutionC_n1_[i])/deltaT; //-solutionC_n_[i](solutionC_n1_[i]-solutionC_n_[i])/deltaT;//
288 }
289
290
291 double time = this->getTimeStep()+deltaT;
292
293 double subIterationTolerance = 1.e-7;
294
295 // immer speicher und wenn es konvergiert, dann zur history machen
296 double historyUpdated[] = {1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 1., 1., 1.512656, 1.512656, 1., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
297 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 1., 1., 1.512656, 1.512656, 1., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
298 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 1., 1., 1.512656, 1.512656, 1., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
299 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 1., 1., 1.512656, 1.512656, 1., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
300
301 AceGenInterface::DeformationDiffusionSmoothMuscleActiveGrowthReorientationTetrahedra3D10 elem(positions, displacements, concentrations, accelerations, rates, &domainData[0], &history[0], subIterationTolerance, deltaT, time, this->iCode_,this->getGlobalElementID());
302 int errorCode = elem.compute();
303
304 TEUCHOS_TEST_FOR_EXCEPTION(errorCode == 2, std::runtime_error, "Subiteration Fail in Element" << this->getGlobalElementID() );
305
306
307 double *residuumRint = elem.getResiduumVectorRint();
308
309 // Write residuumRint to a file named residuumRint.txt
310 // myfile.open ("outputs.txt",ios::app);
311 // myfile << "residuumRint: ";
312 // for(int i=0; i< 30 ; i++){
313 // myfile << residuumRint[i] << " ";
314 // }
315 // myfile << "\n";
316 // myfile.close();
317
318 double *residuumRDyn = elem.getResiduumVectorRdyn();
319
320 for(int i=0; i< 30 ; i++){
321 (*this->rhsVec_)[i] = residuumRint[i]; //+residuumRDyn[i];
322 }
323
324 double *residuumRc = elem.getResiduumVectorRc();
325
326 for(int i=0; i< 10 ; i++){
327 (*this->rhsVec_)[i+30] = residuumRc[i];
328
329 }
330
331
332
333#endif
334
335}
336
353
354
355template <class SC,class LO, class GO, class NO>
356void AssembleFE_SCI_SMC_Active_Growth_Reorientation<SC,LO,GO,NO>::assemble_SCI_SMC_Active_Growth_Reorientation(SmallMatrixPtr_Type &elementMatrix){
357 // double deltat=this->getTimeIncrement();
358 // std::vector<double> deltat(1);
359 // deltat[0]=this->getTimeIncrement();
360#ifdef FEDD_HAVE_ACEGENINTERFACE
361
362 double deltaT=this->getTimeIncrement();
363
364 double positions[30];
365 int count = 0;
366 for(int i=0;i<10;i++){
367 for(int j=0;j<3;j++){
368 positions[count] = this->getNodesRefConfig()[i][j];
369 count++;
370 }
371 }
372 double displacements[30];
373 for(int i = 0; i < 30; i++)
374 {
375 displacements[i]= (*this->solution_)[i];
376
377 }
378 std::vector<double> history(this->historyLength_,0.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, 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
379 //cout << "History_ " ;
380 for(int i = 0; i < this->historyLength_; i++){
381 history[i] = this->history_[i];
382 //cout << history[i] << ", ";
383 }
384 //cout << endl;
385
386
387 double concentrations[10];
388 for(int i = 0; i < 10; i++)
389 {
390 concentrations[i]=(*this->solution_)[i+30];
391 solutionC_n1_[i]=(*this->solution_)[i+30]; // in each newtonstep solution for n+1 is updated.
392 }
393
394 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};
395
396
397 std::vector<double> domainData = {this->fA_, this->lambdaC50_, this->gamma3_, this->lambdaBarCDotMax_, this->lambdaBarCDotMin_, this->gamma2_, this->gamma1_, this->eta_, this->ca50_, this->k2_, this->k5_,
398 this->k3_, this->k4_, this->k7_, this->kappa_, this->beta1_, this->muA_, this->beta2_, this->alpha2_, this->alpha3_, this->alpha1_, this->alpha4_, this->alpha5_,
399 this->gamma6_, this->lambdaP50_, this->kDotMin_, this->zeta1_, this->kDotMax_, this->gamma4_, this->lambdaBarDotPMin_, this->lambdaBarDotPMax_, this->gamma5_, this->zeta2_,
400 this->DeltaLambdaBarPMin_, this->p1_, this->p3_, this->c50_, this->d0_, this->m_, this->activeStartTime_, this->kEtaPlus_, this->mEtaPlus_, this->growthStartTime_,
401 this->reorientationStartTime_, this->growthEndTime_, this->reorientationEndTime_, this->kThetaPlus_, this->kThetaMinus_, this->mThetaPlus_, this->mThetaMinus_, this->thetaPlus1_,
402 this->thetaPlus2_, this->thetaPlus3_, this->thetaMinus1_, this->thetaMinus2_, this->thetaMinus3_, this->kMin_, this->rho_};
403
404
405 double rates[10];
406 for(int i=0; i<10 ; i++){
407 rates[i] =(this->solutionC_n1_[i]-this->solutionC_n_[i]) / deltaT;//
408 }
409 // ##########################
410
411 // history [in] Vector of history variables [Order: LambdaBarC1, LambdaBarC2, nA1, nA2, nB1, nB2, nC1, nC2, nD1, nD2, LambdaA1, LambdaA2] (The length must be equal to number of history variables per gauss point * number of gauss points)
412
413 double time = this->getTimeStep()+deltaT;
414 double subIterationTolerance = 1.e-7;
415
416 // immer speicher und wenn es konvergiert, dann zur history machen
417 // 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
418
419 AceGenInterface::DeformationDiffusionSmoothMuscleActiveGrowthReorientationTetrahedra3D10 elem(positions, displacements, concentrations, accelerations, rates, &domainData[0], &history[0], subIterationTolerance, deltaT, time, this->iCode_,this->getGlobalElementID());
420 int errorCode = elem.compute();
421 TEUCHOS_TEST_FOR_EXCEPTION(errorCode == 2, std::runtime_error, "Subiteration Fail in Element" << this->getGlobalElementID() );
422
423 double *historyUpdated = elem.getHistoryUpdated();
424
425 double** stiffnessMatrixKuu = elem.getStiffnessMatrixKuu();
426 double** stiffnessMatrixKuc = elem.getStiffnessMatrixKuc();
427 double** stiffnessMatrixKcu = elem.getStiffnessMatrixKcu();
428 double** stiffnessMatrixKcc = elem.getStiffnessMatrixKcc();
429 double** massMatrixMc = elem.getMassMatrixMc();
430
431 for(int i=0; i< this->historyLength_; i++){
432 this->historyUpdated_[i] = historyUpdated[i];
433 }
434
435 for(int i=0; i< 30; i++){
436 for(int j=0; j<30; j++){
437 //if(std::fabs(stiffnessMatrixKuu[i][j]) > 1e7)
438 // cout << " !!! Sus entry Kuu [" << i << "][" << j << "] " << stiffnessMatrixKuu[i][j] << endl;
439
440 (*elementMatrix)[i][j]=stiffnessMatrixKuu[i][j];
441 }
442
443 }
444
445 for(int i=0; i< 30; i++){
446 for(int j=0; j<10; j++){
447 //if(std::fabs(stiffnessMatrixKuc[i][j]) > 1e7)
448 // cout << " !!! Sus entry Kuc [" << i << "][" << j << "] " << stiffnessMatrixKuc[i][j] << endl;
449
450 (*elementMatrix)[i][j+30]=stiffnessMatrixKuc[i][j];
451 }
452 }
453 for(int i=0; i< 10; i++){
454 for(int j=0; j<30; j++){
455 //if(std::fabs(stiffnessMatrixKcu[i][j]) > 1e7)
456 // cout << " !!! Sus entry Kcu [" << i << "][" << j << "] " << stiffnessMatrixKcu[i][j] << endl;
457
458 (*elementMatrix)[i+30][j]=stiffnessMatrixKcu[i][j];
459 }
460 }
461 for(int i=0; i< 10; i++){
462 for(int j=0; j<10; j++){
463 //if(std::fabs(massMatrixMc[i][j]) > 1e5 || std::fabs(stiffnessMatrixKcc[i][j]) > 1e5 )
464 // cout << " !!! Sus entry Mass [" << i << "][" << j << "] " << massMatrixMc[i][j] << " or stiff Kcc " << stiffnessMatrixKcc[i][j] << endl;
465
466 (*elementMatrix)[i+30][j+30] =stiffnessMatrixKcc[i][j] +(1./deltaT)*massMatrixMc[i][j]; //
467
468 //(*massMatrix_)[i][j] = massMatrixMc[i][j];
469 }
470 }
471#endif
472
473
474}
475
476} // namespace FEDD
477#endif // AssembleFE_SCI_SMC_Active_Growth_Reorientation_DEF_hpp
FEDD::AssembleFE_SCI_SMC_Active_Growth_Reorientation::assembleJacobian
void assembleJacobian() override
Assemble the element Jacobian matrix.
Definition AssembleFE_SCI_SMC_Active_Growth_Reorientation_def.hpp:192
FEDD::AssembleFE_SCI_SMC_Active_Growth_Reorientation::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_Active_Growth_Reorientation_def.hpp:202
FEDD::AssembleFE_SCI_SMC_Active_Growth_Reorientation::assembleRHS
void assembleRHS() override
Assemble the element right hand side vector.
Definition AssembleFE_SCI_SMC_Active_Growth_Reorientation_def.hpp:226
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