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

Last change on this file since 12706 was 12706, checked in by Mathieu Morlighem, 13 years ago

merged trunk-jpl and trunk for revision 12703
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 "../../../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 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}
Note: See TracBrowser for help on using the repository browser.