source: issm/trunk/src/c/classes/objects/Materials/Matice.cpp@ 13975

Last change on this file since 13975 was 13975, checked in by Mathieu Morlighem, 12 years ago

merged trunk-jpl and trunk for revision 13974
File size: 20.9 KB
Line 
1/*!\file Matice.c
2 * \brief: implementation of the Matice object
3 */
4
5#ifdef HAVE_CONFIG_H
6 #include <config.h>
7#else
8#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
9#endif
10
11#include <stdio.h>
12#include <string.h>
13#include "../../classes.h"
14#include "../../../EnumDefinitions/EnumDefinitions.h"
15#include "../../../shared/shared.h"
16#include "../../../include/include.h"
17
18/*Matice constructors and destructor*/
19/*FUNCTION Matice::Matice(){{{*/
20Matice::Matice(){
21 this->inputs=NULL;
22 this->helement=NULL;
23 return;
24}
25/*}}}*/
26/*FUNCTION Matice::Matice(int id, int index, IoModel* iomodel, int num_vertices){{{*/
27Matice::Matice(int matice_mid,int index, IoModel* iomodel){
28
29 /*Intermediaries:*/
30 int matice_eid;
31
32 /*Initialize id*/
33 this->mid=matice_mid;
34
35 /*Initialize inputs*/
36 this->inputs=new Inputs();
37
38 /*Initialize inputs from IoModel*/
39 this->InputUpdateFromIoModel(index,iomodel);
40
41 /*Hooks: */
42 matice_eid=index+1;
43 this->helement=new Hook(&matice_eid,1);
44
45 return;
46
47}
48/*}}}*/
49/*FUNCTION Matice::~Matice(){{{*/
50Matice::~Matice(){
51 delete helement;
52 delete inputs;
53 return;
54}
55/*}}}*/
56
57/*Object virtual functions definitions:*/
58/*FUNCTION Matice::Echo {{{*/
59void Matice::Echo(void){
60
61 _printLine_("Matice:");
62 _printLine_(" mid: " << mid);
63 _printLine_(" inputs:");
64 inputs->Echo();
65 _printLine_(" element:");
66 helement->Echo();
67}
68/*}}}*/
69/*FUNCTION Matice::DeepEcho {{{*/
70void Matice::DeepEcho(void){
71
72 _printLine_("Matice:");
73 _printLine_(" mid: " << mid);
74 _printLine_(" inputs:");
75 inputs->DeepEcho();
76 _printLine_(" element:");
77 helement->Echo();
78}
79/*}}}*/
80/*FUNCTION Matice::Id {{{*/
81int Matice::Id(void){ return mid; }
82/*}}}*/
83/*FUNCTION Matice::ObjectEnum{{{*/
84int Matice::ObjectEnum(void){
85
86 return MaticeEnum;
87
88}
89/*}}}*/
90/*FUNCTION Matice::copy {{{*/
91Object* Matice::copy() {
92
93 /*Output*/
94 Matice* matice=NULL;
95
96 /*Initialize output*/
97 matice=new Matice();
98
99 /*copy fields: */
100 matice->mid=this->mid;
101 matice->helement=(Hook*)this->helement->copy();
102 if(this->inputs) matice->inputs=(Inputs*)this->inputs->Copy();
103 else matice->inputs=new Inputs();
104
105 return matice;
106}
107/*}}}*/
108
109/*Matice management*/
110/*FUNCTION Matice::Configure {{{*/
111void Matice::Configure(Elements* elementsin){
112
113 /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective
114 * datasets, using internal ids and offsets hidden in hooks: */
115 helement->configure(elementsin);
116}
117/*}}}*/
118/*FUNCTION Matice::SetCurrentConfiguration {{{*/
119void Matice::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){
120
121}
122/*}}}*/
123/*FUNCTION Matice::GetA {{{*/
124IssmDouble Matice::GetA(){
125 /*
126 * A = 1/B^n
127 */
128
129 IssmDouble B,n;
130
131 inputs->GetInputAverage(&B,MaterialsRheologyBEnum);
132 n=this->GetN();
133
134 return pow(B,-n);
135}
136/*}}}*/
137/*FUNCTION Matice::GetB {{{*/
138IssmDouble Matice::GetB(){
139
140 /*Output*/
141 IssmDouble B;
142
143 inputs->GetInputAverage(&B,MaterialsRheologyBEnum);
144 return B;
145}
146/*}}}*/
147/*FUNCTION Matice::GetBbar {{{*/
148IssmDouble Matice::GetBbar(){
149
150 /*Output*/
151 IssmDouble Bbar;
152
153 inputs->GetInputAverage(&Bbar,MaterialsRheologyBbarEnum);
154 return Bbar;
155}
156/*}}}*/
157/*FUNCTION Matice::GetN {{{*/
158IssmDouble Matice::GetN(){
159
160 /*Output*/
161 IssmDouble n;
162
163 inputs->GetInputAverage(&n,MaterialsRheologyNEnum);
164 return n;
165}
166/*}}}*/
167/*FUNCTION Matice::GetVectorFromInputs{{{*/
168void Matice::GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum){
169
170 /*Intermediaries*/
171 Element *element= NULL;
172
173 /*Recover element*/
174 element=(Element*)helement->delivers();
175
176 /*Check that input_enum is a material input*/
177 if (!IsInput(input_enum)) return;
178
179 switch(element->ObjectEnum()){
180
181 case TriaEnum:{
182
183 /*Prepare index list*/
184 int doflist1[3];
185 for(int i=0;i<3;i++) doflist1[i]=((Tria*)element)->nodes[i]->GetVertexPid();
186
187 /*Get input (either in element or material)*/
188 Input* input=inputs->GetInput(input_enum);
189 if(!input) _error_("Input " << EnumToStringx(input_enum) << " not found in material");
190
191 /*We found the enum. Use its values to fill into the vector, using the vertices ids: */
192 input->GetVectorFromInputs(vector,&doflist1[0]);}
193 break;
194
195 default: _error_("element " << EnumToStringx(element->ObjectEnum()) << " not implemented yet");
196 }
197}
198/*}}}*/
199/*FUNCTION Matice::GetViscosity2d {{{*/
200void Matice::GetViscosity2d(IssmDouble* pviscosity, IssmDouble* epsilon){
201 /*From a string tensor and a material object, return viscosity, using Glen's flow law.
202 B
203 viscosity= -------------------------------------------------------------------
204 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]
205
206 where viscosity is the viscotiy, B the flow law parameter , (u,v) the velocity
207 vector, and n the flow law exponent.
208
209 If epsilon is NULL, it means this is the first time SystemMatrices is being run, and we
210 return 10^14, initial viscosity.
211 */
212
213 /*output: */
214 IssmDouble viscosity;
215
216 /*input strain rate: */
217 IssmDouble exx,eyy,exy;
218
219 /*Intermediary: */
220 IssmDouble A,e;
221 IssmDouble B,n;
222
223 /*Get B and n*/
224 B=GetBbar();
225 n=GetN();
226
227 if (n==1){
228 /*Viscous behaviour! viscosity=B: */
229 viscosity=B/2;
230 }
231 else{
232 if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0)){
233 viscosity=0.5*pow((IssmDouble)10,(IssmDouble)14);
234 }
235 else{
236 /*Retrive strain rate components: */
237 exx=*(epsilon+0);
238 eyy=*(epsilon+1);
239 exy=*(epsilon+2);
240
241 /*Build viscosity: viscosity=B/(2*A^e) */
242 A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+exx*eyy;
243 if(A==0){
244 /*Maxiviscositym viscosity for 0 shear areas: */
245 viscosity=2.5*pow(10.,17.);
246 }
247 else{
248 e=(n-1)/(2*n);
249 viscosity=B/(2*pow(A,e));
250 }
251 }
252 }
253
254 /*Checks in debugging mode*/
255 if(viscosity<=0) _error_("Negative viscosity");
256 _assert_(B>0);
257 _assert_(n>0);
258
259 /*Return: */
260 *pviscosity=viscosity;
261}
262/*}}}*/
263/*FUNCTION Matice::GetViscosity3d {{{*/
264void Matice::GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* epsilon){
265
266 /*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]:
267 *
268 * B
269 * viscosity3d= -------------------------------------------------------------------
270 * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]
271 *
272 * where mu is the viscotiy, B the flow law parameter , (u,v) the velocity
273 * vector, and n the flow law exponent.
274 *
275 * If epsilon is NULL, it means this is the first time Emg is being run, and we
276 * return g, initial viscosity.
277 */
278
279 /*output: */
280 IssmDouble viscosity3d;
281
282 /*input strain rate: */
283 IssmDouble exx,eyy,exy,exz,eyz;
284
285 /*Intermediaries: */
286 IssmDouble A,e;
287 IssmDouble B,n;
288
289 /*Get B and n*/
290 B=GetB();
291 n=GetN();
292
293 if (n==1){
294 /*Viscous behaviour! viscosity3d=B: */
295 viscosity3d=B/2;
296 }
297 else{
298 if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) &&
299 (epsilon[3]==0) && (epsilon[4]==0)){
300 viscosity3d=0.5*pow((IssmDouble)10,(IssmDouble)14);
301 }
302 else{
303
304 /*Retrive strain rate components: */
305 exx=*(epsilon+0);
306 eyy=*(epsilon+1);
307 exy=*(epsilon+2);
308 exz=*(epsilon+3);
309 eyz=*(epsilon+4);
310
311 /*Build viscosity: viscosity3d=2*B/(2*A^e) */
312 A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+pow(exz,2)+pow(eyz,2)+exx*eyy;
313 if(A==0){
314 /*Maxiviscosity3dm viscosity for 0 shear areas: */
315 viscosity3d=2.25*pow((IssmDouble)10,(IssmDouble)17);
316 }
317 else{
318 e=(n-1)/2/n;
319
320 viscosity3d=B/(2*pow(A,e));
321 }
322 }
323 }
324
325 /*Checks in debugging mode*/
326 if(viscosity3d<=0) _error_("Negative viscosity");
327 _assert_(B>0);
328 _assert_(n>0);
329
330 /*Assign output pointers:*/
331 *pviscosity3d=viscosity3d;
332}
333/*}}}*/
334/*FUNCTION Matice::GetViscosity3dStokes {{{*/
335void Matice::GetViscosity3dStokes(IssmDouble* pviscosity3d, IssmDouble* epsilon){
336 /*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]:
337 *
338 * B
339 * viscosity3d= -------------------------------------------------------------------
340 * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]
341 *
342 * where mu is the viscotiy, B the flow law parameter , (u,v) the velocity
343 * vector, and n the flow law exponent.
344 *
345 * If epsilon is NULL, it means this is the first time Emg is being run, and we
346 * return g, initial viscosity.
347 */
348
349 /*output: */
350 IssmDouble viscosity3d;
351
352 /*input strain rate: */
353 IssmDouble exx,eyy,exy,exz,eyz,ezz;
354
355 /*Intermediaries: */
356 IssmDouble A,e;
357 IssmDouble B,n;
358 IssmDouble eps0;
359
360 /*Get B and n*/
361 eps0=pow((IssmDouble)10,(IssmDouble)-27);
362 B=GetB();
363 n=GetN();
364
365 if (n==1){
366 /*Viscous behaviour! viscosity3d=B: */
367 viscosity3d=B/2;
368 }
369 else{
370 if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) &&
371 (epsilon[3]==0) && (epsilon[4]==0) && (epsilon[5]==0)){
372 viscosity3d=0.5*pow((IssmDouble)10,(IssmDouble)14);
373 }
374 else{
375
376 /*Retrive strain rate components: */
377 exx=*(epsilon+0);
378 eyy=*(epsilon+1);
379 ezz=*(epsilon+2); //not used
380 exy=*(epsilon+3);
381 exz=*(epsilon+4);
382 eyz=*(epsilon+5);
383
384 /*Build viscosity: viscosity3d=B/(2*A^e) */
385 A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+pow(exz,2)+pow(eyz,2)+exx*eyy+pow(eps0,2);
386 if(A==0){
387 /*Maxiviscosity3dm viscosity for 0 shear areas: */
388 viscosity3d=2.25*pow((IssmDouble)10,(IssmDouble)17);
389 }
390 else{
391 e=(n-1)/2/n;
392 viscosity3d=B/(2*pow(A,e));
393 }
394 }
395 }
396
397 /*Checks in debugging mode*/
398 if(viscosity3d<=0) _error_("Negative viscosity");
399 _assert_(B>0);
400 _assert_(n>0);
401
402 /*Assign output pointers:*/
403 *pviscosity3d=viscosity3d;
404}
405/*}}}*/
406/*FUNCTION Matice::GetViscosityComplement {{{*/
407void Matice::GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){
408 /*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]:
409 *
410 * 1
411 * viscosity= -------------------------------------------------------------------
412 * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]
413 *
414 * If epsilon is NULL, it means this is the first time Gradjb is being run, and we
415 * return mu20, initial viscosity.
416 */
417
418 /*output: */
419 IssmDouble viscosity_complement;
420
421 /*input strain rate: */
422 IssmDouble exx,eyy,exy;
423
424 /*Intermediary value A and exponent e: */
425 IssmDouble A,e;
426 IssmDouble B,n;
427
428 /*Get B and n*/
429 B=GetBbar();
430 n=GetN();
431
432 if(epsilon){
433 exx=*(epsilon+0);
434 eyy=*(epsilon+1);
435 exy=*(epsilon+2);
436
437 /*Build viscosity: mu2=B/(2*A^e) */
438 A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+exx*eyy;
439 if(A==0){
440 /*Maximum viscosity_complement for 0 shear areas: */
441 viscosity_complement=2.25*pow((IssmDouble)10,(IssmDouble)17);
442 }
443 else{
444 e=(n-1)/(2*n);
445
446 viscosity_complement=1/(2*pow(A,e));
447 }
448 }
449 else{
450 viscosity_complement=4.5*pow((IssmDouble)10,(IssmDouble)17);
451 }
452
453 /*Checks in debugging mode*/
454 _assert_(B>0);
455 _assert_(n>0);
456 _assert_(viscosity_complement>0);
457
458 /*Return: */
459 *pviscosity_complement=viscosity_complement;
460}
461/*}}}*/
462/*FUNCTION Matice::GetViscosityDerivativeEpsSquare{{{*/
463void Matice::GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* epsilon){
464
465 /*output: */
466 IssmDouble mu_prime;
467 IssmDouble mu,n,eff2;
468
469 /*input strain rate: */
470 IssmDouble exx,eyy,exy,exz,eyz;
471
472 /*Get visocisty and n*/
473 GetViscosity3d(&mu,epsilon);
474 n=GetN();
475
476 if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) &&
477 (epsilon[3]==0) && (epsilon[4]==0)){
478 mu_prime=0.5*pow((IssmDouble)10,(IssmDouble)14);
479 }
480 else{
481 /*Retrive strain rate components: */
482 exx=epsilon[0];
483 eyy=epsilon[1];
484 exy=epsilon[2];
485 exz=epsilon[3];
486 eyz=epsilon[4];
487 eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy + exz*exz + eyz*eyz;
488
489 mu_prime=(1-n)/(2*n) * mu/eff2;
490 }
491
492 /*Assign output pointers:*/
493 *pmu_prime=mu_prime;
494}
495/*}}}*/
496/*FUNCTION Matice::GetViscosity2dDerivativeEpsSquare{{{*/
497void Matice::GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* epsilon){
498
499 /*output: */
500 IssmDouble mu_prime;
501 IssmDouble mu,n,eff2;
502
503 /*input strain rate: */
504 IssmDouble exx,eyy,exy;
505
506 /*Get visocisty and n*/
507 GetViscosity2d(&mu,epsilon);
508 n=GetN();
509
510 if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0)){
511 mu_prime=0.5*pow((IssmDouble)10,(IssmDouble)14);
512 }
513 else{
514 /*Retrive strain rate components: */
515 exx=epsilon[0];
516 eyy=epsilon[1];
517 exy=epsilon[2];
518 eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy ;
519
520 mu_prime=(1-n)/(2*n) * mu/eff2;
521 }
522
523 /*Assign output pointers:*/
524 *pmu_prime=mu_prime;
525}
526/*}}}*/
527/*FUNCTION Matice::InputDuplicate{{{*/
528void Matice::InputDuplicate(int original_enum,int new_enum){
529
530 /*Call inputs method*/
531 if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum);
532
533}
534/*}}}*/
535/*FUNCTION Matice::InputUpdateFromVector(IssmDouble* vector, int name, int type) {{{*/
536void Matice::InputUpdateFromVector(IssmDouble* vector, int name, int type){
537
538 /*Intermediaries*/
539 Element *element = NULL;
540
541 /*Recover element*/
542 element=(Element*)helement->delivers();
543
544 /*Check that name is an element input*/
545 if (!IsInput(name)) return;
546
547 switch(type){
548
549 case VertexEnum:
550
551 switch(element->ObjectEnum()){
552
553 case TriaEnum: {
554 IssmDouble values[3];
555 for (int i=0;i<3;i++) values[i]=vector[((Tria*)element)->nodes[i]->GetVertexPid()];
556 this->inputs->AddInput(new TriaP1Input(name,values));
557 return;
558 }
559 default: _error_("element " << EnumToStringx(element->ObjectEnum()) << " not implemented yet");
560 }
561 default: _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
562 }
563}
564/*}}}*/
565/*FUNCTION Matice::InputUpdateFromVector(int* vector, int name, int type) {{{*/
566void Matice::InputUpdateFromVector(int* vector, int name, int type){
567 /*Nothing updated yet*/
568}
569/*}}}*/
570/*FUNCTION Matice::InputUpdateFromVector(bool* vector, int name, int type) {{{*/
571void Matice::InputUpdateFromVector(bool* vector, int name, int type){
572 /*Nothing updated yet*/
573}
574/*}}}*/
575/*FUNCTION Matice::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type) {{{*/
576void Matice::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){
577
578 /*Intermediaries*/
579 Element *element = NULL;
580 Parameters* parameters= NULL;
581 int dim;
582
583 /*Recover element*/
584 element=(Element*)helement->delivers();
585
586 /*Check that name is an element input*/
587 if (!IsInput(name)) return;
588
589 switch(type){
590
591 case VertexEnum:
592
593 switch(element->ObjectEnum()){
594
595 case TriaEnum: {
596 IssmDouble values[3];
597 for (int i=0;i<3;i++) values[i]=vector[((Tria*)element)->nodes[i]->GetVertexSid()]; //index into serial oriented vector
598 this->inputs->AddInput(new TriaP1Input(name,values));
599 /*Special case for rheology B in 2D: Pourave land for this solution{{{*/
600 if(name==MaterialsRheologyBEnum){
601 /*Are we in 2D?:*/
602 if(element->ObjectEnum()==TriaEnum){
603 parameters=((Tria*)(element))->parameters;
604 }
605 else{
606 parameters=((Penta*)(element))->parameters;
607 }
608 parameters->FindParam(&dim,MeshDimensionEnum);
609 if(dim==2){
610 /*Dupliacte rheology input: */
611 this->inputs->AddInput(new TriaP1Input(MaterialsRheologyBbarEnum,values));
612 }
613 }
614 /*}}}*/
615 return;
616 }
617 default: _error_("element " << EnumToStringx(element->ObjectEnum()) << " not implemented yet");
618 }
619 default: _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
620 }
621
622}
623/*}}}*/
624/*FUNCTION Matice::InputUpdateFromMatrixDakota(int* vector, int name, int type) {{{*/
625void Matice::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols,int name, int type){
626 /*Nothing updated yet*/
627}
628/*}}}*/
629/*FUNCTION Matice::InputUpdateFromVectorDakota(int* vector, int name, int type) {{{*/
630void Matice::InputUpdateFromVectorDakota(int* vector, int name, int type){
631 /*Nothing updated yet*/
632}
633/*}}}*/
634/*FUNCTION Matice::InputUpdateFromVectorDakota(bool* vector, int name, int type) {{{*/
635void Matice::InputUpdateFromVectorDakota(bool* vector, int name, int type){
636 /*Nothing updated yet*/
637}
638/*}}}*/
639/*FUNCTION Matice::InputUpdateFromConstant(IssmDouble constant, int name) {{{*/
640void Matice::InputUpdateFromConstant(IssmDouble constant, int name){
641 /*Nothing updated yet*/
642}
643/*}}}*/
644/*FUNCTION Matice::InputUpdateFromConstant(int constant, int name) {{{*/
645void Matice::InputUpdateFromConstant(int constant, int name){
646 /*Nothing updated yet*/
647}
648/*}}}*/
649/*FUNCTION Matice::InputUpdateFromConstant(bool constant, int name) {{{*/
650void Matice::InputUpdateFromConstant(bool constant, int name){
651 /*Nothing updated yet*/
652}
653/*}}}*/
654/*FUNCTION Matice::InputUpdateFromSolution{{{*/
655void Matice::InputUpdateFromSolution(IssmDouble* solution){
656 /*Nothing updated yet*/
657}
658/*}}}*/
659/*FUNCTION Matice::InputUpdateFromIoModel{{{*/
660void Matice::InputUpdateFromIoModel(int index, IoModel* iomodel){
661
662 int i,j;
663
664 int dim;
665 bool control_analysis;
666 int num_control_type;
667
668 /*Fetch parameters: */
669 iomodel->Constant(&dim,MeshDimensionEnum);
670 iomodel->Constant(&control_analysis,InversionIscontrolEnum);
671 if(control_analysis) iomodel->Constant(&num_control_type,InversionNumControlParametersEnum);
672
673 /*if 2d*/
674 if(dim==2){
675
676 /*Intermediaries*/
677 const int num_vertices = 3; //Tria has 3 vertices
678 IssmDouble nodeinputs[num_vertices];
679 IssmDouble cmmininputs[num_vertices];
680 IssmDouble cmmaxinputs[num_vertices];
681
682 /*Get B*/
683 if (iomodel->Data(MaterialsRheologyBEnum)) {
684 for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyBEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)];
685 this->inputs->AddInput(new TriaP1Input(MaterialsRheologyBbarEnum,nodeinputs));
686 }
687
688 /*Get n*/
689 if (iomodel->Data(MaterialsRheologyNEnum)) {
690 for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyNEnum)[index];
691 this->inputs->AddInput(new TriaP1Input(MaterialsRheologyNEnum,nodeinputs));
692 }
693
694 /*Control Inputs*/
695 #ifdef _HAVE_CONTROL_
696 if (control_analysis && iomodel->Data(InversionControlParametersEnum)){
697 for(i=0;i<num_control_type;i++){
698 switch(reCast<int>(iomodel->Data(InversionControlParametersEnum)[i])){
699 case MaterialsRheologyBbarEnum:
700 if (iomodel->Data(MaterialsRheologyBEnum)){
701 _assert_(iomodel->Data(MaterialsRheologyBEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum));
702 for(j=0;j<num_vertices;j++)nodeinputs[j]=iomodel->Data(MaterialsRheologyBEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)];
703 for(j=0;j<num_vertices;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
704 for(j=0;j<num_vertices;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
705 this->inputs->AddInput(new ControlInput(MaterialsRheologyBbarEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
706 }
707 break;
708 }
709 }
710 }
711 #endif
712 }
713
714 /*if 3d*/
715 #ifdef _HAVE_3D_
716 else if(dim==3){
717
718 /*Intermediaries*/
719 const int num_vertices = 6; //Penta has 6 vertices
720 IssmDouble nodeinputs[num_vertices];
721 IssmDouble cmmininputs[num_vertices];
722 IssmDouble cmmaxinputs[num_vertices];
723
724 /*Get B*/
725 if (iomodel->Data(MaterialsRheologyBEnum)) {
726 for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyBEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)];
727 this->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,nodeinputs));
728 }
729
730 /*Get n*/
731 if (iomodel->Data(MaterialsRheologyNEnum)) {
732 for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyNEnum)[index];
733 this->inputs->AddInput(new PentaP1Input(MaterialsRheologyNEnum,nodeinputs));
734 }
735
736 /*Control Inputs*/
737 #ifdef _HAVE_CONTROL_
738 if (control_analysis && iomodel->Data(InversionControlParametersEnum)){
739 for(i=0;i<num_control_type;i++){
740 switch(reCast<int>(iomodel->Data(InversionControlParametersEnum)[i])){
741 case MaterialsRheologyBbarEnum:
742 if (iomodel->Data(MaterialsRheologyBEnum)){
743 _assert_(iomodel->Data(MaterialsRheologyBEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum));
744 for(j=0;j<num_vertices;j++)nodeinputs[j]=iomodel->Data(MaterialsRheologyBEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)];
745 for(j=0;j<num_vertices;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
746 for(j=0;j<num_vertices;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
747 this->inputs->AddInput(new ControlInput(MaterialsRheologyBEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
748 }
749 break;
750 }
751 }
752 }
753 #endif
754 }
755 #endif
756 else{
757 _error_("Mesh type not supported yet!");
758 }
759
760 return;
761}
762/*}}}*/
763/*FUNCTION Matice::IsInput{{{*/
764bool Matice::IsInput(int name){
765 if (
766 name==MaterialsRheologyBEnum ||
767 name==MaterialsRheologyBbarEnum ||
768 name==MaterialsRheologyNEnum
769 ){
770 return true;
771 }
772 else return false;
773}
774/*}}}*/
Note: See TracBrowser for help on using the repository browser.