source: issm/oecreview/Archive/15392-16133/ISSM-15483-15484.diff@ 16134

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

Added Archive/15392-16133
File size: 12.6 KB
RevLine 
[16134]1Index: ../trunk-jpl/src/c/modules/ModelProcessorx/Prognostic/CreateConstraintsPrognostic.cpp
2===================================================================
3--- ../trunk-jpl/src/c/modules/ModelProcessorx/Prognostic/CreateConstraintsPrognostic.cpp (revision 15483)
4+++ ../trunk-jpl/src/c/modules/ModelProcessorx/Prognostic/CreateConstraintsPrognostic.cpp (revision 15484)
5@@ -15,7 +15,7 @@
6
7 /*Do not add constraints in DG, they are weakly imposed*/
8 if(stabilization!=3){
9- IoModelToConstraintsx(constraints,iomodel,PrognosticSpcthicknessEnum,PrognosticAnalysisEnum);
10+ IoModelToConstraintsx(constraints,iomodel,PrognosticSpcthicknessEnum,PrognosticAnalysisEnum,P1Enum);
11 }
12
13 /*Assign output pointer: */
14Index: ../trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CreateConstraintsHydrologyShreve.cpp
15===================================================================
16--- ../trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CreateConstraintsHydrologyShreve.cpp (revision 15483)
17+++ ../trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CreateConstraintsHydrologyShreve.cpp (revision 15484)
18@@ -19,7 +19,7 @@
19
20 if(hydrology_model!=HydrologyshreveEnum) return;
21
22- IoModelToConstraintsx(constraints,iomodel,HydrologyshreveSpcwatercolumnEnum,HydrologyShreveAnalysisEnum);
23+ IoModelToConstraintsx(constraints,iomodel,HydrologyshreveSpcwatercolumnEnum,HydrologyShreveAnalysisEnum,P1Enum);
24
25 /*Assign output pointer: */
26 *pconstraints=constraints;
27Index: ../trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp
28===================================================================
29--- ../trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (revision 15483)
30+++ ../trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (revision 15484)
31@@ -10,18 +10,18 @@
32 void CreateConstraintsDiagnosticHoriz(Constraints** pconstraints, IoModel* iomodel){
33
34 /*Intermediary*/
35- int i,j;
36- int count;
37- IssmDouble yts;
38- IssmDouble g;
39- IssmDouble rho_ice;
40- IssmDouble stokesreconditioning;
41- bool isstokes,isl1l2,ismacayealpattyn;
42- int fe_ssa;
43- bool spcpresent=false;
44- int Mx,Nx;
45- int My,Ny;
46- int Mz,Nz;
47+ int i,j;
48+ int count;
49+ IssmDouble yts;
50+ IssmDouble g;
51+ IssmDouble rho_ice;
52+ IssmDouble stokesreconditioning;
53+ bool isstokes,isl1l2,ismacayealpattyn;
54+ int fe_ssa;
55+ bool spcpresent = false;
56+ int Mx,Nx;
57+ int My,Ny;
58+ int Mz,Nz;
59 IssmDouble *spcvx = NULL;
60 IssmDouble *spcvy = NULL;
61 IssmDouble *spcvz = NULL;
62@@ -328,18 +328,18 @@
63
64 if(!xIsNan<IssmDouble>(spcvx[v1]) && !xIsNan<IssmDouble>(spcvx[v2])){
65 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,
66- 1,(spcvx[v1]+spcvx[v2])/yts,DiagnosticHorizAnalysisEnum));
67+ 1,(spcvx[v1]+spcvx[v2])/(2.*yts),DiagnosticHorizAnalysisEnum));
68 count++;
69 }
70 if(!xIsNan<IssmDouble>(spcvy[v1]) && !xIsNan<IssmDouble>(spcvy[v2])){
71 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,
72- 2,(spcvy[v1]+spcvy[v2])/yts,DiagnosticHorizAnalysisEnum));
73+ 2,(spcvy[v1]+spcvy[v2])/(2.*yts),DiagnosticHorizAnalysisEnum));
74 count++;
75 }
76 if (reCast<int,IssmDouble>(vertices_type[v1])==StokesApproximationEnum || (reCast<int,IssmDouble>(vertices_type[v1])==NoneApproximationEnum)){
77 if(!xIsNan<IssmDouble>(spcvz[v1]) && !xIsNan<IssmDouble>(spcvz[v2])){
78 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,
79- 3,(spcvz[v1]+spcvz[v2])/yts,DiagnosticHorizAnalysisEnum));
80+ 3,(spcvz[v1]+spcvz[v2])/(2.*yts),DiagnosticHorizAnalysisEnum));
81 count++;
82 }
83 }
84Index: ../trunk-jpl/src/c/modules/ModelProcessorx/Balancethickness/CreateConstraintsBalancethickness.cpp
85===================================================================
86--- ../trunk-jpl/src/c/modules/ModelProcessorx/Balancethickness/CreateConstraintsBalancethickness.cpp (revision 15483)
87+++ ../trunk-jpl/src/c/modules/ModelProcessorx/Balancethickness/CreateConstraintsBalancethickness.cpp (revision 15484)
88@@ -15,7 +15,7 @@
89
90 /*Do not add constraints in DG*/
91 if(stabilization!=3){
92- IoModelToConstraintsx(constraints,iomodel,BalancethicknessSpcthicknessEnum,BalancethicknessAnalysisEnum);
93+ IoModelToConstraintsx(constraints,iomodel,BalancethicknessSpcthicknessEnum,BalancethicknessAnalysisEnum,P1Enum);
94 }
95
96 /*Assign output pointer: */
97Index: ../trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CreateConstraintsHydrologyDCEfficient.cpp
98===================================================================
99--- ../trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CreateConstraintsHydrologyDCEfficient.cpp (revision 15483)
100+++ ../trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CreateConstraintsHydrologyDCEfficient.cpp (revision 15484)
101@@ -23,7 +23,7 @@
102 iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum);
103 if(!isefficientlayer) return;
104
105- IoModelToConstraintsx(constraints,iomodel,HydrologydcSpceplHeadEnum,HydrologyDCEfficientAnalysisEnum);
106+ IoModelToConstraintsx(constraints,iomodel,HydrologydcSpceplHeadEnum,HydrologyDCEfficientAnalysisEnum,P1Enum);
107
108 /*Assign output pointer: */
109 *pconstraints=constraints;
110Index: ../trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateConstraintsHydrologyDCInefficient.cpp
111===================================================================
112--- ../trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateConstraintsHydrologyDCInefficient.cpp (revision 15483)
113+++ ../trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateConstraintsHydrologyDCInefficient.cpp (revision 15484)
114@@ -18,7 +18,7 @@
115 iomodel->Constant(&hydrology_model,HydrologyModelEnum);
116 if(hydrology_model!=HydrologydcEnum) return;
117
118- IoModelToConstraintsx(constraints,iomodel,HydrologydcSpcsedimentHeadEnum,HydrologyDCInefficientAnalysisEnum);
119+ IoModelToConstraintsx(constraints,iomodel,HydrologydcSpcsedimentHeadEnum,HydrologyDCInefficientAnalysisEnum,P1Enum);
120
121 /*Assign output pointer: */
122 *pconstraints=constraints;
123Index: ../trunk-jpl/src/c/modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp
124===================================================================
125--- ../trunk-jpl/src/c/modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp (revision 15483)
126+++ ../trunk-jpl/src/c/modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp (revision 15484)
127@@ -15,7 +15,7 @@
128
129 /*Only 3d mesh supported*/
130 if(iomodel->dim==3){
131- IoModelToConstraintsx(constraints,iomodel,ThermalSpctemperatureEnum,ThermalAnalysisEnum);
132+ IoModelToConstraintsx(constraints,iomodel,ThermalSpctemperatureEnum,ThermalAnalysisEnum,P1Enum);
133 }
134
135 /*Assign output pointer: */
136Index: ../trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.h
137===================================================================
138--- ../trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.h (revision 15483)
139+++ ../trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.h (revision 15484)
140@@ -7,6 +7,6 @@
141 #include "../../classes/classes.h"
142
143 /* local prototypes: */
144-void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,int vector_enum,int analysis_type);
145+void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,int vector_enum,int analysis_type,int finite_element);
146
147 #endif /* _IOMODELTOELEMENTINPUTX_H */
148Index: ../trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp
149===================================================================
150--- ../trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp (revision 15483)
151+++ ../trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp (revision 15484)
152@@ -5,55 +5,67 @@
153 #include "./IoModelToConstraintsx.h"
154 #include "../../shared/shared.h"
155 #include "../../toolkits/toolkits.h"
156+#include "../ModelProcessorx/ModelProcessorx.h"
157
158-void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,int vector_enum,int analysis_type){
159+void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,int vector_enum,int analysis_type,int finite_element){
160
161 /*intermediary: */
162- int i,j;
163- IssmDouble yts;
164- bool transient = false;
165+ int i,j,count;
166+ IssmDouble yts;
167 FILE *fid = NULL;
168- int code = 0;
169- int vector_layout = 0;
170+ int code,vector_layout;
171 IssmDouble *times = NULL;
172 IssmDouble *values = NULL;
173 bool spcpresent = false;
174- int count = 0;
175
176+ /*P2 finite elements*/
177+ int v1,v2;
178+ bool *my_edges = NULL;
179+
180 /*variables being fetched: */
181- IssmDouble *IssmDoublevector = NULL;
182- int M,N;
183+ IssmDouble *spcdata = NULL;
184+ int M,N;
185
186 /*Fetch parameters: */
187 iomodel->Constant(&yts,ConstantsYtsEnum);
188
189 /*First of, find the record for the enum, and get code of data type: */
190 fid=iomodel->SetFilePointerToData(&code, &vector_layout,vector_enum);
191-
192 if(code!=7)_error_("expecting a IssmDouble vector for constraints with enum " << EnumToStringx(vector_enum));
193 if(vector_layout!=1)_error_("expecting a nodal vector for constraints with enum " << EnumToStringx(vector_enum));
194
195 /*Fetch vector:*/
196- iomodel->FetchData(&IssmDoublevector,&M,&N,vector_enum);
197+ iomodel->FetchData(&spcdata,&M,&N,vector_enum);
198
199- /*Transient or static?:*/
200+ /*Partition edges if we are using P2 finite elements*/
201+ if(finite_element==P2Enum){
202+ EdgesPartitioning(&my_edges,iomodel);
203+ }
204+
205 if(M==iomodel->numberofvertices){
206- /*static: just create Constraints objects*/
207+ /*Static constraint*/
208 count=0;
209-
210- /*Create Constraints from x,y,z: */
211 for (i=0;i<iomodel->numberofvertices;i++){
212-
213- /*keep only this partition's nodes:*/
214 if((iomodel->my_vertices[i])){
215-
216- if (!xIsNan<IssmDouble>(IssmDoublevector[i])){
217-
218- constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,IssmDoublevector[i],analysis_type));
219+ if (!xIsNan<IssmDouble>(spcdata[i])){
220+ constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcdata[i],analysis_type));
221 count++;
222 }
223 }
224 }
225+ if(finite_element==P2Enum){
226+ for(i=0;i<iomodel->numberofedges;i++){
227+ if(my_edges[i]){
228+ v1 = iomodel->edges[4*i+0]-1;
229+ v2 = iomodel->edges[4*i+1]-1;
230+ if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
231+ constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,
232+ 1,(spcdata[v1]+spcdata[v2])/2.,analysis_type));
233+ count++;
234+ }
235+ }
236+ }
237+ }
238 }
239 else if (M==(iomodel->numberofvertices+1)){
240 /*transient: create transient SpcTransient objects. Same logic, except we need to retrieve
241@@ -62,21 +74,17 @@
242
243 /*figure out times: */
244 times=xNew<IssmDouble>(N);
245- for(j=0;j<N;j++){
246- times[j]=IssmDoublevector[(M-1)*N+j]*yts;
247- }
248+ for(j=0;j<N;j++) times[j]=spcdata[(M-1)*N+j]*yts;
249
250- /*Create constraints from x,y,z: */
251+ /*Create constraints: */
252 for (i=0;i<iomodel->numberofvertices;i++){
253-
254- /*keep only this partition's nodes:*/
255 if((iomodel->my_vertices[i])){
256
257 /*figure out times and values: */
258 values=xNew<IssmDouble>(N);
259 spcpresent=false;
260 for(j=0;j<N;j++){
261- values[j]=IssmDoublevector[i*N+j];
262+ values[j]=spcdata[i*N+j];
263 if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
264 }
265
266@@ -87,13 +95,32 @@
267 xDelete<IssmDouble>(values);
268 }
269 }
270+ if(finite_element==P2Enum){
271+ for(i=0;i<iomodel->numberofedges;i++){
272+ if(my_edges[i]){
273+ v1 = iomodel->edges[4*i+0]-1;
274+ v2 = iomodel->edges[4*i+1]-1;
275+ values=xNew<IssmDouble>(N);
276+ spcpresent=false;
277+ for(j=0;j<N;j++){
278+ values[j]=(spcdata[v1*N+j]+spcdata[v2*N+j])/2.;
279+ if(!xIsNan<IssmDouble>(values[j])) spcpresent=true; //NaN means no spc by default
280+ }
281+ if(spcpresent){
282+ constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,1,
283+ N,times,values,analysis_type));
284+ count++;
285+ }
286+ }
287+ }
288+ }
289 }
290 else{
291 _error_("Size of field " << EnumToStringx(vector_enum) << " not supported");
292 }
293
294 /*Free ressources:*/
295- xDelete<IssmDouble>(IssmDoublevector);
296+ xDelete<IssmDouble>(spcdata);
297 xDelete<IssmDouble>(times);
298 xDelete<IssmDouble>(values);
299 }
Note: See TracBrowser for help on using the repository browser.