Index: ../trunk-jpl/src/c/modules/ModelProcessorx/Prognostic/CreateConstraintsPrognostic.cpp =================================================================== --- ../trunk-jpl/src/c/modules/ModelProcessorx/Prognostic/CreateConstraintsPrognostic.cpp (revision 15483) +++ ../trunk-jpl/src/c/modules/ModelProcessorx/Prognostic/CreateConstraintsPrognostic.cpp (revision 15484) @@ -15,7 +15,7 @@ /*Do not add constraints in DG, they are weakly imposed*/ if(stabilization!=3){ - IoModelToConstraintsx(constraints,iomodel,PrognosticSpcthicknessEnum,PrognosticAnalysisEnum); + IoModelToConstraintsx(constraints,iomodel,PrognosticSpcthicknessEnum,PrognosticAnalysisEnum,P1Enum); } /*Assign output pointer: */ Index: ../trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CreateConstraintsHydrologyShreve.cpp =================================================================== --- ../trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CreateConstraintsHydrologyShreve.cpp (revision 15483) +++ ../trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CreateConstraintsHydrologyShreve.cpp (revision 15484) @@ -19,7 +19,7 @@ if(hydrology_model!=HydrologyshreveEnum) return; - IoModelToConstraintsx(constraints,iomodel,HydrologyshreveSpcwatercolumnEnum,HydrologyShreveAnalysisEnum); + IoModelToConstraintsx(constraints,iomodel,HydrologyshreveSpcwatercolumnEnum,HydrologyShreveAnalysisEnum,P1Enum); /*Assign output pointer: */ *pconstraints=constraints; Index: ../trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp =================================================================== --- ../trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (revision 15483) +++ ../trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (revision 15484) @@ -10,18 +10,18 @@ void CreateConstraintsDiagnosticHoriz(Constraints** pconstraints, IoModel* iomodel){ /*Intermediary*/ - int i,j; - int count; - IssmDouble yts; - IssmDouble g; - IssmDouble rho_ice; - IssmDouble stokesreconditioning; - bool isstokes,isl1l2,ismacayealpattyn; - int fe_ssa; - bool spcpresent=false; - int Mx,Nx; - int My,Ny; - int Mz,Nz; + int i,j; + int count; + IssmDouble yts; + IssmDouble g; + IssmDouble rho_ice; + IssmDouble stokesreconditioning; + bool isstokes,isl1l2,ismacayealpattyn; + int fe_ssa; + bool spcpresent = false; + int Mx,Nx; + int My,Ny; + int Mz,Nz; IssmDouble *spcvx = NULL; IssmDouble *spcvy = NULL; IssmDouble *spcvz = NULL; @@ -328,18 +328,18 @@ if(!xIsNan(spcvx[v1]) && !xIsNan(spcvx[v2])){ constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1, - 1,(spcvx[v1]+spcvx[v2])/yts,DiagnosticHorizAnalysisEnum)); + 1,(spcvx[v1]+spcvx[v2])/(2.*yts),DiagnosticHorizAnalysisEnum)); count++; } if(!xIsNan(spcvy[v1]) && !xIsNan(spcvy[v2])){ constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1, - 2,(spcvy[v1]+spcvy[v2])/yts,DiagnosticHorizAnalysisEnum)); + 2,(spcvy[v1]+spcvy[v2])/(2.*yts),DiagnosticHorizAnalysisEnum)); count++; } if (reCast(vertices_type[v1])==StokesApproximationEnum || (reCast(vertices_type[v1])==NoneApproximationEnum)){ if(!xIsNan(spcvz[v1]) && !xIsNan(spcvz[v2])){ constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1, - 3,(spcvz[v1]+spcvz[v2])/yts,DiagnosticHorizAnalysisEnum)); + 3,(spcvz[v1]+spcvz[v2])/(2.*yts),DiagnosticHorizAnalysisEnum)); count++; } } Index: ../trunk-jpl/src/c/modules/ModelProcessorx/Balancethickness/CreateConstraintsBalancethickness.cpp =================================================================== --- ../trunk-jpl/src/c/modules/ModelProcessorx/Balancethickness/CreateConstraintsBalancethickness.cpp (revision 15483) +++ ../trunk-jpl/src/c/modules/ModelProcessorx/Balancethickness/CreateConstraintsBalancethickness.cpp (revision 15484) @@ -15,7 +15,7 @@ /*Do not add constraints in DG*/ if(stabilization!=3){ - IoModelToConstraintsx(constraints,iomodel,BalancethicknessSpcthicknessEnum,BalancethicknessAnalysisEnum); + IoModelToConstraintsx(constraints,iomodel,BalancethicknessSpcthicknessEnum,BalancethicknessAnalysisEnum,P1Enum); } /*Assign output pointer: */ Index: ../trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CreateConstraintsHydrologyDCEfficient.cpp =================================================================== --- ../trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CreateConstraintsHydrologyDCEfficient.cpp (revision 15483) +++ ../trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CreateConstraintsHydrologyDCEfficient.cpp (revision 15484) @@ -23,7 +23,7 @@ iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum); if(!isefficientlayer) return; - IoModelToConstraintsx(constraints,iomodel,HydrologydcSpceplHeadEnum,HydrologyDCEfficientAnalysisEnum); + IoModelToConstraintsx(constraints,iomodel,HydrologydcSpceplHeadEnum,HydrologyDCEfficientAnalysisEnum,P1Enum); /*Assign output pointer: */ *pconstraints=constraints; Index: ../trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateConstraintsHydrologyDCInefficient.cpp =================================================================== --- ../trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateConstraintsHydrologyDCInefficient.cpp (revision 15483) +++ ../trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateConstraintsHydrologyDCInefficient.cpp (revision 15484) @@ -18,7 +18,7 @@ iomodel->Constant(&hydrology_model,HydrologyModelEnum); if(hydrology_model!=HydrologydcEnum) return; - IoModelToConstraintsx(constraints,iomodel,HydrologydcSpcsedimentHeadEnum,HydrologyDCInefficientAnalysisEnum); + IoModelToConstraintsx(constraints,iomodel,HydrologydcSpcsedimentHeadEnum,HydrologyDCInefficientAnalysisEnum,P1Enum); /*Assign output pointer: */ *pconstraints=constraints; Index: ../trunk-jpl/src/c/modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp =================================================================== --- ../trunk-jpl/src/c/modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp (revision 15483) +++ ../trunk-jpl/src/c/modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp (revision 15484) @@ -15,7 +15,7 @@ /*Only 3d mesh supported*/ if(iomodel->dim==3){ - IoModelToConstraintsx(constraints,iomodel,ThermalSpctemperatureEnum,ThermalAnalysisEnum); + IoModelToConstraintsx(constraints,iomodel,ThermalSpctemperatureEnum,ThermalAnalysisEnum,P1Enum); } /*Assign output pointer: */ Index: ../trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.h =================================================================== --- ../trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.h (revision 15483) +++ ../trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.h (revision 15484) @@ -7,6 +7,6 @@ #include "../../classes/classes.h" /* local prototypes: */ -void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,int vector_enum,int analysis_type); +void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,int vector_enum,int analysis_type,int finite_element); #endif /* _IOMODELTOELEMENTINPUTX_H */ Index: ../trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp =================================================================== --- ../trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp (revision 15483) +++ ../trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp (revision 15484) @@ -5,55 +5,67 @@ #include "./IoModelToConstraintsx.h" #include "../../shared/shared.h" #include "../../toolkits/toolkits.h" +#include "../ModelProcessorx/ModelProcessorx.h" -void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,int vector_enum,int analysis_type){ +void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,int vector_enum,int analysis_type,int finite_element){ /*intermediary: */ - int i,j; - IssmDouble yts; - bool transient = false; + int i,j,count; + IssmDouble yts; FILE *fid = NULL; - int code = 0; - int vector_layout = 0; + int code,vector_layout; IssmDouble *times = NULL; IssmDouble *values = NULL; bool spcpresent = false; - int count = 0; + /*P2 finite elements*/ + int v1,v2; + bool *my_edges = NULL; + /*variables being fetched: */ - IssmDouble *IssmDoublevector = NULL; - int M,N; + IssmDouble *spcdata = NULL; + int M,N; /*Fetch parameters: */ iomodel->Constant(&yts,ConstantsYtsEnum); /*First of, find the record for the enum, and get code of data type: */ fid=iomodel->SetFilePointerToData(&code, &vector_layout,vector_enum); - if(code!=7)_error_("expecting a IssmDouble vector for constraints with enum " << EnumToStringx(vector_enum)); if(vector_layout!=1)_error_("expecting a nodal vector for constraints with enum " << EnumToStringx(vector_enum)); /*Fetch vector:*/ - iomodel->FetchData(&IssmDoublevector,&M,&N,vector_enum); + iomodel->FetchData(&spcdata,&M,&N,vector_enum); - /*Transient or static?:*/ + /*Partition edges if we are using P2 finite elements*/ + if(finite_element==P2Enum){ + EdgesPartitioning(&my_edges,iomodel); + } + if(M==iomodel->numberofvertices){ - /*static: just create Constraints objects*/ + /*Static constraint*/ count=0; - - /*Create Constraints from x,y,z: */ for (i=0;inumberofvertices;i++){ - - /*keep only this partition's nodes:*/ if((iomodel->my_vertices[i])){ - - if (!xIsNan(IssmDoublevector[i])){ - - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,IssmDoublevector[i],analysis_type)); + if (!xIsNan(spcdata[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcdata[i],analysis_type)); count++; } } } + if(finite_element==P2Enum){ + for(i=0;inumberofedges;i++){ + if(my_edges[i]){ + v1 = iomodel->edges[4*i+0]-1; + v2 = iomodel->edges[4*i+1]-1; + if(!xIsNan(spcdata[v1]) && !xIsNan(spcdata[v2])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1, + 1,(spcdata[v1]+spcdata[v2])/2.,analysis_type)); + count++; + } + } + } + } } else if (M==(iomodel->numberofvertices+1)){ /*transient: create transient SpcTransient objects. Same logic, except we need to retrieve @@ -62,21 +74,17 @@ /*figure out times: */ times=xNew(N); - for(j=0;jnumberofvertices;i++){ - - /*keep only this partition's nodes:*/ if((iomodel->my_vertices[i])){ /*figure out times and values: */ values=xNew(N); spcpresent=false; for(j=0;j(values[j]))spcpresent=true; //NaN means no spc by default } @@ -87,13 +95,32 @@ xDelete(values); } } + if(finite_element==P2Enum){ + for(i=0;inumberofedges;i++){ + if(my_edges[i]){ + v1 = iomodel->edges[4*i+0]-1; + v2 = iomodel->edges[4*i+1]-1; + values=xNew(N); + spcpresent=false; + for(j=0;j(values[j])) spcpresent=true; //NaN means no spc by default + } + if(spcpresent){ + constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,1, + N,times,values,analysis_type)); + count++; + } + } + } + } } else{ _error_("Size of field " << EnumToStringx(vector_enum) << " not supported"); } /*Free ressources:*/ - xDelete(IssmDoublevector); + xDelete(spcdata); xDelete(times); xDelete(values); }