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