1 | /*
|
---|
2 | * CreateConstraintsEnthalpy.c:
|
---|
3 | */
|
---|
4 |
|
---|
5 | #include "../../../Container/Container.h"
|
---|
6 | #include "../../../io/io.h"
|
---|
7 | #include "../../../toolkits/toolkits.h"
|
---|
8 | #include "../../../EnumDefinitions/EnumDefinitions.h"
|
---|
9 | #include "../../../objects/objects.h"
|
---|
10 | #include "../../../shared/shared.h"
|
---|
11 | #include "../../../include/include.h"
|
---|
12 | #include "../ModelProcessorx.h"
|
---|
13 |
|
---|
14 | void CreateConstraintsEnthalpy(Constraints** pconstraints, IoModel* iomodel){
|
---|
15 |
|
---|
16 | /*Intermediary*/
|
---|
17 | int i,j;
|
---|
18 | int count;
|
---|
19 | int dim;
|
---|
20 | int M,N;
|
---|
21 | int numberofvertices;
|
---|
22 | bool spcpresent=false;
|
---|
23 | double heatcapacity;
|
---|
24 | double referencetemperature;
|
---|
25 |
|
---|
26 | /*Output*/
|
---|
27 | IssmDouble *spcvector = NULL;
|
---|
28 | IssmDouble* times=NULL;
|
---|
29 | IssmDouble* values=NULL;
|
---|
30 | Constraints* constraints = NULL;
|
---|
31 |
|
---|
32 | /*Fetch parameters: */
|
---|
33 | iomodel->Constant(&dim,MeshDimensionEnum);
|
---|
34 | iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum);
|
---|
35 | iomodel->Constant(&heatcapacity,MaterialsHeatcapacityEnum);
|
---|
36 | iomodel->Constant(&referencetemperature,ConstantsReferencetemperatureEnum);
|
---|
37 |
|
---|
38 | /*Recover pointer: */
|
---|
39 | constraints=*pconstraints;
|
---|
40 |
|
---|
41 | /*Create constraints if they do not exist yet*/
|
---|
42 | if(!constraints) constraints = new Constraints();
|
---|
43 |
|
---|
44 | /*return if 2d mesh*/
|
---|
45 | if (dim==2){
|
---|
46 | *pconstraints=constraints;
|
---|
47 | return;
|
---|
48 | }
|
---|
49 |
|
---|
50 | /*Fetch data: */
|
---|
51 | iomodel->FetchData(&spcvector,&M,&N,ThermalSpctemperatureEnum);
|
---|
52 |
|
---|
53 | //FIX ME: SHOULD USE IOMODELCREATECONSTRAINTS
|
---|
54 | /*Transient or static?:*/
|
---|
55 | if(M==numberofvertices){
|
---|
56 | /*static: just create Constraints objects*/
|
---|
57 | count=0;
|
---|
58 |
|
---|
59 | for (i=0;i<numberofvertices;i++){
|
---|
60 | /*keep only this partition's nodes:*/
|
---|
61 | if((iomodel->my_vertices[i])){
|
---|
62 |
|
---|
63 | if (!xIsNan<IssmDouble>(spcvector[i])){
|
---|
64 |
|
---|
65 | constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,heatcapacity*(spcvector[i]-referencetemperature),EnthalpyAnalysisEnum));
|
---|
66 | count++;
|
---|
67 |
|
---|
68 | }
|
---|
69 | }
|
---|
70 | }
|
---|
71 | }
|
---|
72 | else if (M==(numberofvertices+1)){
|
---|
73 | /*transient: create transient SpcTransient objects. Same logic, except we need to retrieve
|
---|
74 | * various times and values to initialize an SpcTransient object: */
|
---|
75 | count=0;
|
---|
76 |
|
---|
77 | /*figure out times: */
|
---|
78 | times=xNew<IssmDouble>(N);
|
---|
79 | for(j=0;j<N;j++){
|
---|
80 | times[j]=spcvector[(M-1)*N+j];
|
---|
81 | }
|
---|
82 | /*unit conversion: */
|
---|
83 | UnitConversion(times,N,ExtToIuEnum,TimeEnum);
|
---|
84 |
|
---|
85 | /*Create constraints from x,y,z: */
|
---|
86 | for (i=0;i<numberofvertices;i++){
|
---|
87 |
|
---|
88 | /*keep only this partition's nodes:*/
|
---|
89 | if((iomodel->my_vertices[i])){
|
---|
90 |
|
---|
91 | /*figure out times and values: */
|
---|
92 | values=xNew<IssmDouble>(N);
|
---|
93 | spcpresent=false;
|
---|
94 | for(j=0;j<N;j++){
|
---|
95 | values[j]=heatcapacity*(spcvector[i*N+j]-referencetemperature);
|
---|
96 | if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
|
---|
97 | }
|
---|
98 |
|
---|
99 | if(spcpresent){
|
---|
100 | constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,N,times,values,EnthalpyAnalysisEnum));
|
---|
101 | count++;
|
---|
102 | }
|
---|
103 | xDelete<IssmDouble>(values);
|
---|
104 | }
|
---|
105 | }
|
---|
106 | }
|
---|
107 | else{
|
---|
108 | _error2_("Size of field " << EnumToStringx(ThermalSpctemperatureEnum) << " not supported");
|
---|
109 | }
|
---|
110 |
|
---|
111 | /*Free ressources:*/
|
---|
112 | xDelete<IssmDouble>(spcvector);
|
---|
113 | xDelete<IssmDouble>(times);
|
---|
114 | xDelete<IssmDouble>(values);
|
---|
115 |
|
---|
116 | /*Assign output pointer: */
|
---|
117 | *pconstraints=constraints;
|
---|
118 | }
|
---|