source: issm/trunk/src/c/modules/ModelProcessorx/Enthalpy/CreateConstraintsEnthalpy.cpp@ 15396

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

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