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