Index: ../trunk-jpl/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp =================================================================== --- ../trunk-jpl/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp (revision 23045) +++ ../trunk-jpl/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp (revision 23046) @@ -1,6 +1,6 @@ /*!\file: InterpFromGridToMeshx.cpp * \brief "c" core code for interpolating values from a structured grid. - */ + */ #ifdef HAVE_CONFIG_H #include @@ -24,9 +24,6 @@ double* y=NULL; int i; - // TEST - _printf_("\r N: "<x1 && y2>y1); _assert_(x<=x2 && x>=x1 && y<=y2 && y>=y1); - return + return +Q11*(x2-x)*(y2-y)/((x2-x1)*(y2-y1)) +Q21*(x-x1)*(y2-y)/((x2-x1)*(y2-y1)) +Q12*(x2-x)*(y-y1)/((x2-x1)*(y2-y1)) @@ -300,7 +297,7 @@ * | | | * | | | * y1 x--------x---------x Q21 - * x1 xm x2 + * x1 xm x2 * */ /*Checks*/ Index: ../trunk-jpl/src/wrappers/InterpFromGridToMesh/InterpFromGridToMesh.js =================================================================== --- ../trunk-jpl/src/wrappers/InterpFromGridToMesh/InterpFromGridToMesh.js (revision 23045) +++ ../trunk-jpl/src/wrappers/InterpFromGridToMesh/InterpFromGridToMesh.js (revision 23046) @@ -5,7 +5,7 @@ * * Usage: * var data_mesh=InterpFromGridToMesh(xIn,yIn,dataIn,xMeshIn,yMeshIn,defaultValue,interpType);\ - * + * * xIn,yIn : coordinates of matrix data. (x and y must be in increasing order) * dataIn : matrix holding the data to be interpolated onto the mesh * xMeshIn,yMeshIn : coordinates of the points onto which we interpolate @@ -50,13 +50,13 @@ var y = {}; var yMesh = {}; //}}} - + /* Dynamic allocations */ //{{{ - + /* Input */ @@ -67,7 +67,7 @@ dxHeap = new Uint8Array(Module.HEAPU8.buffer, dxPtr, nx); dxHeap.set(new Uint8Array(dx.buffer)); x = dxHeap.byteOffset; - + dy = new Float64Array(yIn); ny = dy.length * dy.BYTES_PER_ELEMENT; dyPtr = Module._malloc(ny); @@ -74,7 +74,7 @@ dyHeap = new Uint8Array(Module.HEAPU8.buffer, dyPtr, ny); dyHeap.set(new Uint8Array(dy.buffer)); y = dyHeap.byteOffset; - + ddata = new Float64Array(dataIn); ndata = ddata.length * ddata.BYTES_PER_ELEMENT; ddataPtr = Module._malloc(ndata); @@ -81,7 +81,7 @@ ddataHeap = new Uint8Array(Module.HEAPU8.buffer, ddataPtr, ndata); ddataHeap.set(new Uint8Array(ddata.buffer)); data = ddataHeap.byteOffset; - + dxMesh = new Float64Array(xMeshIn); nxMesh = dxMesh.length * dxMesh.BYTES_PER_ELEMENT; dxMeshPtr = Module._malloc(nxMesh); @@ -88,7 +88,7 @@ dxMeshHeap = new Uint8Array(Module.HEAPU8.buffer, dxMeshPtr, nxMesh); dxMeshHeap.set(new Uint8Array(dxMesh.buffer)); xMesh = dxMeshHeap.byteOffset; - + dyMesh = new Float64Array(yMeshIn); nyMesh = dyMesh.length * dyMesh.BYTES_PER_ELEMENT; dyMeshPtr = Module._malloc(nyMesh); @@ -95,12 +95,11 @@ dyMeshHeap = new Uint8Array(Module.HEAPU8.buffer, dyMeshPtr, nyMesh); dyMeshHeap.set(new Uint8Array(dyMesh.buffer)); yMesh = dyMeshHeap.byteOffset; - + nods = xMeshIn.length; meshNumRows = xMeshIn.length; - console.log('InterpFromGridToMesh.js: meshNumRows => ' + meshNumRows); - - + + /* Retrieve interpolation type */ @@ -111,23 +110,23 @@ interpType = 'bilinear'; } //}}} - + /* Output */ pdataMesh = Module._malloc(4); //}}} - + //}}} - - + + /* Declare InterpFromGridToMesh module */ //{{{ InterpFromGridToMeshModule = Module.cwrap( - 'InterpFromGridToMeshModule', - 'number', + 'InterpFromGridToMeshModule', + 'number', [ 'number', // output : pdataMesh 'number', // input : x @@ -144,20 +143,20 @@ ] ); //}}} - - + + /* Call InterpFromGridToMesh module */ //{{{ InterpFromGridToMeshModule( - pdataMesh, - x, - y, - data, - xMesh, - yMesh, - defaultValue, + pdataMesh, + x, + y, + data, + xMesh, + yMesh, + defaultValue, nods, dataNumRowsIn, dataNumColsIn, @@ -165,8 +164,8 @@ interpType ); //}}} - - + + /* Dynamic copying from heap */ @@ -174,8 +173,8 @@ dataMeshPtr = Module.getValue(pdataMesh, 'i32'); dataMesh = Module.HEAPF64.slice(dataMeshPtr / 8, dataMeshPtr / 8 + nods); //}}} - - + + /* Free resources */ @@ -182,7 +181,7 @@ //{{{ Module._free(pdataMesh); //}}} - - + + return dataMesh; } Index: ../trunk-jpl/jenkins/macosx_pine-island =================================================================== --- ../trunk-jpl/jenkins/macosx_pine-island (revision 23045) +++ ../trunk-jpl/jenkins/macosx_pine-island (revision 23046) @@ -4,9 +4,9 @@ #-------------------------------# #MATLAB path -MATLAB_PATH="/Applications/MATLAB_R2017b.app/" +MATLAB_PATH="/Applications/MATLAB_R2015b.app/" -#ISSM CONFIGURATION +#ISSM CONFIGURATION ISSM_CONFIG='--prefix=$ISSM_DIR \ --with-matlab-dir=$MATLAB_PATH \ --with-triangle-dir=$ISSM_DIR/externalpackages/triangle/install \ @@ -53,6 +53,6 @@ #as follows: runme($MATLAB_NROPTIONS). The options must be understandable #by Matlab and runme.m #ex: "'id',[101 102 103]" -## FS +## FS PYTHON_NROPTIONS="" MATLAB_NROPTIONS="'exclude',[701,702,703,435,IdFromString('Dakota')]" Index: ../trunk-jpl/jenkins/macosx_pine-island_static =================================================================== --- ../trunk-jpl/jenkins/macosx_pine-island_static (revision 23045) +++ ../trunk-jpl/jenkins/macosx_pine-island_static (revision 23046) @@ -4,9 +4,9 @@ #-------------------------------# #MATLAB path -MATLAB_PATH="/Applications/MATLAB_R2017b.app/" +MATLAB_PATH="/Applications/MATLAB_R2015b.app/" -#ISSM CONFIGURATION +#ISSM CONFIGURATION ISSM_CONFIG='--prefix=$ISSM_DIR \ --disable-static \ --enable-standalone-executables \ @@ -22,7 +22,7 @@ --with-metis-dir=$ISSM_DIR/externalpackages/petsc/install \ --with-m1qn3-dir=$ISSM_DIR/externalpackages/m1qn3/install \ --with-math77-dir=$ISSM_DIR/externalpackages/math77/install \ - --with-fortran-lib="/usr/local/gfortran/lib/libgfortran.a /usr/local/gfortran/lib/libquadmath.a /usr/local/gfortran/lib/gcc/x86_64-apple-darwin15/6.1.0/libgcc.a" \ + --with-fortran-lib="/usr/local/gfortran/lib/libgfortran.a /usr/local/gfortran/lib/libquadmath.a /usr/local/gfortran/lib/gcc/x86_64-apple-darwin14/5.2.0/libgcc.a" \ --with-numthreads=4' #PYTHON and MATLAB testing @@ -60,6 +60,6 @@ #as follows: runme($MATLAB_NROPTIONS). The options must be understandable #by Matlab and runme.m #ex: "'id',[101 102 103]" -## bamg mesh FS +## bamg mesh FS #PYTHON_NROPTIONS="" #MATLAB_NROPTIONS="'exclude',[119,243,514,701,702,703,435,IdFromString('Dakota')]" Index: ../trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 23045) +++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 23046) @@ -1463,7 +1463,7 @@ /*Figure out maximum across the cluster: */ ISSM_MPI_Reduce(&maxvy,&node_maxvy,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&node_maxvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&node_maxvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); maxvy=node_maxvy; /*Assign output pointers:*/ @@ -1487,7 +1487,7 @@ /*Figure out maximum across the cluster: */ ISSM_MPI_Reduce(&maxvz,&node_maxvz,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&node_maxvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&node_maxvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); maxvz=node_maxvz; /*Assign output pointers:*/ @@ -1511,7 +1511,7 @@ /*Figure out minimum across the cluster: */ ISSM_MPI_Reduce(&minvel,&node_minvel,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&node_minvel,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&node_minvel,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); minvel=node_minvel; /*Assign output pointers:*/ @@ -1535,7 +1535,7 @@ /*Figure out minimum across the cluster: */ ISSM_MPI_Reduce(&minvx,&node_minvx,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&node_minvx,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&node_minvx,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); minvx=node_minvx; /*Assign output pointers:*/ @@ -1559,7 +1559,7 @@ /*Figure out minimum across the cluster: */ ISSM_MPI_Reduce(&minvy,&node_minvy,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&node_minvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&node_minvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); minvy=node_minvy; /*Assign output pointers:*/ @@ -1583,7 +1583,7 @@ /*Figure out minimum across the cluster: */ ISSM_MPI_Reduce(&minvz,&node_minvz,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&node_minvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&node_minvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); minvz=node_minvz; /*Assign output pointers:*/ @@ -1628,7 +1628,7 @@ weights_input->GetInputValue(&weight,gauss,OmegaAbsGradientEnum); omega_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss); - /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ + /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ //J+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight; J+=weight*1/2*pow(dp[0]*dp[0]+dp[1]*dp[1],2)*Jdet*gauss->weight; } @@ -1964,7 +1964,7 @@ else{ IssmDouble* values = xNewZeroInit(nlines*ncols); IssmDouble* allvalues = xNew(nlines*ncols); - + /*Fill-in matrix*/ for(int j=0;jSize();j++){ Element* element=xDynamicCast(this->elements->GetObjectByOffset(j)); @@ -1973,7 +1973,7 @@ /*Gather from all cpus*/ ISSM_MPI_Allreduce((void*)values,(void*)allvalues,ncols*nlines,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); xDelete(values); - + if(save_results)results->AddResult(new GenericExternalResult(results->Size()+1,output_enum,allvalues,nlines,ncols,step,time)); xDelete(allvalues); } @@ -2074,13 +2074,13 @@ case MaterialsRheologyBbarEnum: this->ElementResponsex(responses,MaterialsRheologyBbarEnum); break; case VelEnum: this->ElementResponsex(responses,VelEnum); break; case FrictionCoefficientEnum: NodalValuex(responses, FrictionCoefficientEnum,elements,nodes, vertices, loads, materials, parameters); break; - default: + default: if(response_descriptor_enum>=Outputdefinition1Enum && response_descriptor_enum <=Outputdefinition100Enum){ IssmDouble double_result = OutputDefinitionsResponsex(this,response_descriptor_enum); *responses=double_result; } - else _error_("response descriptor \"" << EnumToStringx(response_descriptor_enum) << "\" not supported yet!"); - break; + else _error_("response descriptor \"" << EnumToStringx(response_descriptor_enum) << "\" not supported yet!"); + break; } } @@ -2211,7 +2211,7 @@ weights_input->GetInputValue(&weight,gauss,ThicknessAbsGradientEnum); thickness_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss); - /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ + /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ J+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight; } @@ -2402,9 +2402,9 @@ Analysis* analysis= EnumToAnalysis(analysis_type); analysis->UpdateConstraints(this); delete analysis; - + /*Second, constraints might be time dependent: */ - SpcNodesx(nodes,constraints,parameters,analysis_type); + SpcNodesx(nodes,constraints,parameters,analysis_type); /*Now, update degrees of freedoms: */ NodesDofx(nodes,parameters,analysis_type); @@ -2415,7 +2415,7 @@ IssmDouble *surface = NULL; IssmDouble *bed = NULL; - + if(VerboseSolution()) _printf0_(" updating vertices positions\n"); /*get vertex vectors for bed and thickness: */ @@ -2454,7 +2454,7 @@ /*}}}*/ /*AMR*/ -#if !defined(_HAVE_ADOLC_) +#if !defined(_HAVE_ADOLC_) void FemModel::ReMesh(void){/*{{{*/ /*Intermediaries*/ @@ -2464,13 +2464,13 @@ int *newelementslist = NULL; int newnumberofvertices = -1; int newnumberofelements = -1; - bool* my_elements = NULL; + bool* my_elements = NULL; int* my_vertices = NULL; int elementswidth = this->GetElementsWidth();//just tria elements in this version int amrtype,basalforcing_model; bool isgroundingline; - /*Branch to specific amr depending on requested method*/ + /*Branch to specific amr depending on requested method*/ this->parameters->FindParam(&amrtype,AmrTypeEnum); switch(amrtype){ #if defined(_HAVE_NEOPZ_) && !defined(_HAVE_ADOLC_) @@ -2483,7 +2483,7 @@ default: _error_("not implemented yet"); } - + /*Partitioning the new mesh. Maybe ElementsAndVerticesPartitioning.cpp could be modified to set this without iomodel.*/ this->ElementsAndVerticesPartitioning(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,&my_elements,&my_vertices); @@ -2514,7 +2514,7 @@ int analysis_enum = this->analysis_type_list[i]; /*As the domain is 2D, it is not necessary to create nodes for this analysis*/ - if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue; + if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue; this->CreateNodes(newnumberofvertices,my_vertices,nodecounter,analysis_enum,new_nodes); this->CreateConstraints(new_vertices,nodecounter,constraintcounter,analysis_enum,new_constraints); @@ -2544,7 +2544,7 @@ analysis_type=this->analysis_type_list[i]; //SetCurrentConfiguration(analysis_type); - this->analysis_counter=i; + this->analysis_counter=i; /*Now, plug analysis_counter and analysis_type inside the parameters: */ this->parameters->SetParam(this->analysis_counter,AnalysisCounterEnum); this->parameters->SetParam(analysis_type,AnalysisTypeEnum); @@ -2561,7 +2561,7 @@ } ConfigureObjectsx(new_elements,this->loads,new_nodes,new_vertices,new_materials,this->parameters); - if(i==0){ + if(i==0){ VerticesDofx(new_vertices,this->parameters); //only call once, we only have one set of vertices } SpcNodesx(new_nodes,new_constraints,this->parameters,analysis_type); @@ -2605,7 +2605,7 @@ /*Insert bedrock from mismip+ setup*/ /*This was used to Misomip project/simulations*/ - + if(VerboseSolution())_printf0_(" call Mismip bedrock adjust module\n"); IssmDouble x,y,bx,by; @@ -2622,7 +2622,7 @@ bx = -150.-728.8*pow(x/300000.,2)+343.91*pow(x/300000.,4)-50.57*pow(x/300000.,6); by = 500./(1.+exp((-2./4000.)*(y-80000./2.-24000.)))+500./(1.+exp((2./4000.)*(y-80000./2.+24000.))); r[i] = max(bx+by,-720.); - } + } /*insert new bedrock*/ element->AddInput(BedEnum,&r[0],P1Enum); /*Cleanup*/ @@ -2635,7 +2635,7 @@ void FemModel::AdjustBaseThicknessAndMask(void){/*{{{*/ if(VerboseSolution())_printf0_(" call adjust base and thickness module\n"); - + int numvertices = this->GetElementsWidth(); IssmDouble rho_water,rho_ice,density,base_float; IssmDouble* phi = xNew(numvertices); @@ -2647,7 +2647,7 @@ for(int el=0;elelements->Size();el++){ Element* element=xDynamicCast(this->elements->GetObjectByOffset(el)); - + element->GetInputListOnVertices(&s[0],SurfaceEnum); element->GetInputListOnVertices(&r[0],BedEnum); element->GetInputListOnVertices(&sl[0],SealevelEnum); @@ -2659,16 +2659,16 @@ /*calculate base floatation (which supports given surface*/ base_float = rho_ice*s[i]/(rho_ice-rho_water); if(r[i]>base_float){ - b[i] = r[i]; - } + b[i] = r[i]; + } else { - b[i] = base_float; - } + b[i] = base_float; + } if(fabs(sl[i])>0) _error_("Sea level value "<AddInput(ThicknessEnum,&h[0],P1Enum); element->AddInput(BaseEnum,&b[0],P1Enum); } - + /*Delete*/ xDelete(phi); xDelete(h); @@ -2704,7 +2704,7 @@ int* P1input_interp = NULL; Vector* input_interpolations = NULL; IssmDouble* input_interpolations_serial = NULL; - int* pos = NULL; + int* pos = NULL; IssmDouble value = 0; /*Figure out how many inputs we have and their respective interpolation*/ @@ -2724,7 +2724,7 @@ P0input_enums = xNew(numP0inputs); P0input_interp = xNew(numP0inputs); P1input_enums = xNew(numP1inputs); - P1input_interp = xNew(numP1inputs); + P1input_interp = xNew(numP1inputs); } numP0inputs = 0; numP1inputs = 0; @@ -2756,23 +2756,23 @@ } } } - - /*Get P0 and P1 inputs over the elements*/ + + /*Get P0 and P1 inputs over the elements*/ pos = xNew(elementswidth); vP0inputs= new Vector(numberofelements*numP0inputs); vP1inputs= new Vector(numberofvertices*numP1inputs); for(int i=0;ielements->Size();i++){ Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - + /*Get P0 inputs*/ for(int j=0;j(element->GetInput(P0input_enums[j])); + TriaInput* input=xDynamicCast(element->GetInput(P0input_enums[j])); input->GetInputAverage(&value); pos[0]=element->Sid()*numP0inputs+j; - /*Insert input in the vector*/ + /*Insert input in the vector*/ vP0inputs->SetValues(1,pos,&value,INS_VAL); } - + /*Get P1 inputs*/ for(int j=0;j(element->GetInput(P1input_enums[j])); @@ -2779,8 +2779,8 @@ pos[0]=element->vertices[0]->Sid()*numP1inputs+j; pos[1]=element->vertices[1]->Sid()*numP1inputs+j; pos[2]=element->vertices[2]->Sid()*numP1inputs+j; - /*Insert input in the vector*/ - vP1inputs->SetValues(elementswidth,pos,input->values,INS_VAL); + /*Insert input in the vector*/ + vP1inputs->SetValues(elementswidth,pos,input->values,INS_VAL); } } @@ -2789,16 +2789,16 @@ vP1inputs->Assemble(); P0inputs=vP0inputs->ToMPISerial(); P1inputs=vP1inputs->ToMPISerial(); - + /*Assign pointers*/ - *pnumP0inputs = numP0inputs; - *pP0inputs = P0inputs; + *pnumP0inputs = numP0inputs; + *pP0inputs = P0inputs; *pP0input_enums = P0input_enums; *pP0input_interp = P0input_interp; - *pnumP1inputs = numP1inputs; - *pP1inputs = P1inputs; + *pnumP1inputs = numP1inputs; + *pP1inputs = P1inputs; *pP1input_enums = P1input_enums; - *pP1input_interp = P1input_interp; + *pP1input_interp = P1input_interp; /*Cleanup*/ delete input_interpolations; @@ -2809,7 +2809,7 @@ } /*}}}*/ void FemModel::InterpolateInputs(Vertices* newfemmodel_vertices,Elements* newfemmodel_elements){/*{{{*/ - + int numberofelements = this->elements->NumberOfElements(); //global, entire old mesh int newnumberofelements = newfemmodel_elements->Size(); //just on the new partition int numberofvertices = this->vertices->NumberOfVertices(); //global, entire old mesh @@ -2825,7 +2825,7 @@ IssmDouble* newP1inputs = NULL; //just on the new partition int* P1input_enums = NULL; int* P1input_interp = NULL; - IssmDouble* values = NULL; + IssmDouble* values = NULL; IssmDouble* vector = NULL; IssmDouble* x = NULL;//global, entire old mesh IssmDouble* y = NULL;//global, entire old mesh @@ -2837,7 +2837,7 @@ IssmDouble* newxc = NULL;//just on the new partition IssmDouble* newyc = NULL;//just on the new partition int* newelementslist = NULL;//just on the new partition - int* sidtoindex = NULL;//global vertices sid to partition index + int* sidtoindex = NULL;//global vertices sid to partition index /*Get old P0 and P1 inputs (entire mesh)*/ this->GetInputs(&numP0inputs,&P0inputs,&P0input_enums,&P0input_interp,&numP1inputs,&P1inputs,&P1input_enums,&P1input_interp); @@ -2870,17 +2870,17 @@ newx,newy,newnumberofvertices,NULL); /*Insert P0 and P1 inputs into the new elements (just on the new partition)*/ - values=xNew(elementswidth); + values=xNew(elementswidth); for(int i=0;iSize();i++){//just on the new partition Element* element=xDynamicCast(newfemmodel_elements->GetObjectByOffset(i)); /*newP0inputs is just on the new partition*/ for(int j=0;jAddInput(new DoubleInput(P0input_enums[j],newP0inputs[i*numP0inputs+j])); break; - case IntInputEnum: + case IntInputEnum: element->AddInput(new IntInput(P0input_enums[j],reCast(newP0inputs[i*numP0inputs+j]))); break; case BoolInputEnum: @@ -2898,7 +2898,7 @@ element->inputs->AddInput(new TriaInput(P1input_enums[j],values,P1Enum)); } } - + /*Cleanup*/ xDelete(P0inputs); xDelete(newP0inputs); @@ -2937,7 +2937,7 @@ int* elementslist = NULL; if(!this->elements || !this->vertices || !this->results || !this->parameters) return; - + parameters->FindParam(&step,StepEnum); parameters->FindParam(&time,TimeEnum); numberofelements=this->elements->NumberOfElements(); @@ -2955,7 +2955,7 @@ this->results->AddResult(new GenericExternalResult(this->results->Size()+1,MeshYEnum, y,numberofvertices,1,step,time)); - + /*Cleanup*/ xDelete(x); xDelete(y); @@ -3016,7 +3016,7 @@ /*Assemble*/ vmasklevelset->Assemble(); - + /*Serialize and set output*/ (*pmasklevelset)=vmasklevelset->ToMPISerial(); @@ -3030,7 +3030,7 @@ void FemModel::CreateVertices(int newnumberofvertices,int newnumberofelements,int elementswidth,int* newelementslist,int* my_vertices,IssmDouble* newx,IssmDouble* newy,IssmDouble* newz,Vertices* vertices){/*{{{*/ /*newelementslist is in Matlab indexing*/ - + /*Creating connectivity table*/ int* connectivity=NULL; connectivity=xNewZeroInit(newnumberofvertices); @@ -3041,12 +3041,12 @@ _assert_(vertexid>0 && vertexid-1id=i+1; newvertex->sid=i; newvertex->pid=UNDEF; @@ -3057,8 +3057,8 @@ newvertex->sigma=0.; newvertex->connectivity=connectivity[i]; newvertex->clone=false;//itapopo check this - vertices->AddObject(newvertex); - } + vertices->AddObject(newvertex); + } } xDelete(connectivity); @@ -3085,13 +3085,13 @@ for(int j=0;jelement_type_list[j]=0; } else newtria->element_type_list=NULL; - + /*Element hook*/ int matpar_id=newnumberofelements+1; //retrieve material parameter id (last pointer in femodel->materials) int material_id=i+1; // retrieve material_id = i+1; /*retrieve vertices ids*/ int* vertex_ids=xNew(elementswidth); - for(int j=0;j(newelementslist[elementswidth*i+j]);//this Hook wants Matlab indexing + for(int j=0;j(newelementslist[elementswidth*i+j]);//this Hook wants Matlab indexing /*Setting the hooks*/ newtria->numanalyses =this->nummodels; newtria->hnodes =new Hook*[this->nummodels]; @@ -3103,8 +3103,8 @@ for(int j=0;jnummodels;j++) newtria->hnodes[j]=NULL; /*Clean up*/ xDelete(vertex_ids); - elements->AddObject(newtria); - } + elements->AddObject(newtria); + } } } @@ -3114,14 +3114,14 @@ /*Just Matice in this version*/ for(int i=0;iAddObject(new Matice(i+1,i,MaticeEnum)); - } + materials->AddObject(new Matice(i+1,i,MaticeEnum)); + } } - + /*Add new constant material property to materials, at the end: */ Matpar *newmatpar=static_cast(this->materials->GetObjectByOffset(this->materials->Size()-1)->copy()); newmatpar->SetMid(newnumberofelements+1); - materials->AddObject(newmatpar);//put it at the end of the materials + materials->AddObject(newmatpar);//put it at the end of the materials } /*}}}*/ void FemModel::CreateNodes(int newnumberofvertices,int* my_vertices,int nodecounter,int analysis_enum,Nodes* nodes){/*{{{*/ @@ -3128,23 +3128,23 @@ int lid=0; for(int j=0;jid=nodecounter+j+1; newnode->sid=j; newnode->lid=lid++; newnode->analysis_enum=analysis_enum; - + /*Initialize coord_system: Identity matrix by default*/ for(int k=0;k<3;k++) for(int l=0;l<3;l++) newnode->coord_system[k][l]=0.0; for(int k=0;k<3;k++) newnode->coord_system[k][k]=1.0; - + /*indexing:*/ newnode->indexingupdate=true; - + Analysis* analysis=EnumToAnalysis(analysis_enum); int *doftypes=NULL; int numdofs=analysis->DofsPerNode(&doftypes,Domain2DhorizontalEnum,SSAApproximationEnum); @@ -3158,7 +3158,7 @@ /*Stressbalance Horiz*/ if(analysis_enum==StressbalanceAnalysisEnum){ - // itapopo this code is rarely used. + // itapopo this code is rarely used. /*Coordinate system provided, convert to coord_system matrix*/ //XZvectorsToCoordinateSystem(&this->coord_system[0][0],&iomodel->Data(StressbalanceReferentialEnum)[j*6]); //_assert_(sqrt( coord_system[0][0]*coord_system[0][0] + coord_system[1][0]*coord_system[1][0]) >1.e-4); @@ -3173,13 +3173,13 @@ if(!femmodel_vertices) _error_("GetMesh: vertices are NULL."); if(!femmodel_elements) _error_("GetMesh: elements are NULL."); - + int numberofvertices = femmodel_vertices->NumberOfVertices(); int numberofelements = femmodel_elements->NumberOfElements(); int elementswidth = this->GetElementsWidth(); // just 2D mesh in this version (just tria elements) IssmDouble* x = NULL; IssmDouble* y = NULL; - IssmDouble* z = NULL; + IssmDouble* z = NULL; int* elementslist = NULL; int* elem_vertices = NULL; IssmDouble *id1 = NULL; @@ -3188,7 +3188,7 @@ /*Get vertices coordinates*/ VertexCoordinatesx(&x,&y,&z,femmodel_vertices,false) ; - + /*Get element vertices*/ elem_vertices = xNew(elementswidth); Vector* vid1= new Vector(numberofelements); @@ -3203,7 +3203,7 @@ vid2->SetValue(element->sid,elem_vertices[1],INS_VAL); vid3->SetValue(element->sid,elem_vertices[2],INS_VAL); } - + /*Assemble*/ vid1->Assemble(); vid2->Assemble(); @@ -3213,7 +3213,7 @@ id1 = vid1->ToMPISerial(); id2 = vid2->ToMPISerial(); id3 = vid3->ToMPISerial(); - + /*Construct elements list*/ elementslist=xNew(numberofelements*elementswidth); if(numberofelements*elementswidth<0) _error_("numberofelements negative."); @@ -3222,7 +3222,7 @@ elementslist[elementswidth*i+1] = reCast(id2[i])+1; //InterpMesh wants Matlab indexing elementslist[elementswidth*i+2] = reCast(id3[i])+1; //InterpMesh wants Matlab indexing } - + /*Assign pointers*/ *px = x; *py = y; @@ -3243,7 +3243,7 @@ if(!femmodel_vertices) _error_("GetMesh: vertices are NULL."); if(!femmodel_elements) _error_("GetMesh: elements are NULL."); - + int numberofvertices = femmodel_vertices->Size(); //number of vertices of this partition int numbertotalofvertices = femmodel_vertices->NumberOfVertices(); //number total of vertices (entire mesh) int numberofelements = femmodel_elements->Size(); //number of elements of this partition @@ -3250,20 +3250,20 @@ int elementswidth = this->GetElementsWidth(); //just 2D mesh in this version (just tria elements) IssmDouble* x = NULL; IssmDouble* y = NULL; - IssmDouble* z = NULL; + IssmDouble* z = NULL; int* elementslist = NULL; int* sidtoindex = NULL; int* elem_vertices = NULL; - + /*Get vertices coordinates of this partition*/ sidtoindex = xNewZeroInit(numbertotalofvertices);//entire mesh, all vertices x = xNew(numberofvertices);//just this partition y = xNew(numberofvertices);//just this partitio; z = xNew(numberofvertices);//just this partitio; - + /*Go through in this partition (vertices)*/ for(int i=0;iGetObjectByOffset(i); + Vertex* vertex=(Vertex*)femmodel_vertices->GetObjectByOffset(i); /*Attention: no spherical coordinates*/ x[i]=vertex->GetX(); y[i]=vertex->GetY(); @@ -3276,7 +3276,7 @@ elem_vertices= xNew(elementswidth); elementslist = xNew(numberofelements*elementswidth); if(numberofelements*elementswidth<0) _error_("numberofelements negative."); - + for(int i=0;i(femmodel_elements->GetObjectByOffset(i)); element->GetVerticesSidList(elem_vertices); @@ -3283,14 +3283,14 @@ elementslist[elementswidth*i+0] = sidtoindex[elem_vertices[0]]+1; //InterpMesh wants Matlab indexing elementslist[elementswidth*i+1] = sidtoindex[elem_vertices[1]]+1; //InterpMesh wants Matlab indexing elementslist[elementswidth*i+2] = sidtoindex[elem_vertices[2]]+1; //InterpMesh wants Matlab indexing - } - + } + /*Assign pointers*/ *px = x; *py = y; *pz = z; *pelementslist = elementslist; //Matlab indexing. InterMesh uses this type. - *psidtoindex = sidtoindex; //it is ncessary to insert inputs + *psidtoindex = sidtoindex; //it is ncessary to insert inputs /*Cleanup*/ xDelete(elem_vertices); @@ -3301,9 +3301,9 @@ /*ATTENTION: JUST SPCVX AND SPCVY*/ /*OTHERS CONSTRAINTS MUST BE IMPLEMENTED*/ if(analysis_enum!=StressbalanceAnalysisEnum) return; - + int numberofnodes_analysistype= this->nodes->NumberOfNodes(analysis_enum); - int dofpernode = 2; //vx and vy + int dofpernode = 2; //vx and vy int numberofcols = dofpernode*2; //to keep dofs and flags in the vspc vector int numberofvertices = this->vertices->NumberOfVertices(); //global, entire old mesh int numberofelements = this->elements->NumberOfElements(); //global, entire old mesh @@ -3327,7 +3327,7 @@ newx=xNew(newnumberofvertices);//just the new partition newy=xNew(newnumberofvertices);//just the new partition for(int i=0;iGetObjectByOffset(i); + Vertex* vertex=(Vertex*)newfemmodel_vertices->GetObjectByOffset(i); /*Attention: no spherical coordinates*/ newx[i]=vertex->GetX(); newy[i]=vertex->GetY(); @@ -3335,7 +3335,7 @@ /*Get spcvx and spcvy of old mesh*/ for(int i=0;iconstraints->Size();i++){ - + Constraint* constraint=(Constraint*)constraints->GetObjectByOffset(i); if(!constraint->InAnalysis(analysis_enum)) _error_("AMR create constraints for "<(constraint); int dof = spcstatic->GetDof(); int node = spcstatic->GetNodeId(); - IssmDouble spcvalue = spcstatic->GetValue(); + IssmDouble spcvalue = spcstatic->GetValue(); int nodeindex = node-1; - + /*vx and vx flag insertion*/ if(dof==0) {//vx vspc->SetValue(nodeindex*numberofcols,spcvalue,INS_VAL); //vx vspc->SetValue(nodeindex*numberofcols+dofpernode,1,INS_VAL);//vxflag - } + } /*vy and vy flag insertion*/ if(dof==1){//vy vspc->SetValue(nodeindex*numberofcols+1,spcvalue,INS_VAL); //vy @@ -3365,7 +3365,7 @@ InterpFromMeshToMesh2dx(&newspc,elementslist,x,y,numberofvertices,numberofelements, spc,numberofvertices,numberofcols, newx,newy,newnumberofvertices,NULL); - + /*Now, insert the interpolated constraints in the data set (constraints)*/ count=0; for(int i=0;iGetObjectByOffset(i); /*spcvy*/ if(!xIsNan(newspc[i*numberofcols+1]) && newspc[i*numberofcols+dofpernode+1]>(1-eps) ){ - newfemmodel_constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+vertex->Sid()+1,1,newspc[i*numberofcols+1],analysis_enum)); + newfemmodel_constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+vertex->Sid()+1,1,newspc[i*numberofcols+1],analysis_enum)); //add count'th spc, on node i+1, setting dof 1 to vx. count++; } @@ -3439,13 +3439,13 @@ int numprocs = IssmComm::GetSize(); bool *my_elements = NULL; int *my_vertices = NULL; - - _assert_(newnumberofvertices>0); - _assert_(newnumberofelements>0); + + _assert_(newnumberofvertices>0); + _assert_(newnumberofelements>0); epart=xNew(newnumberofelements); npart=xNew(newnumberofvertices); index=xNew(elementswidth*newnumberofelements); - + for (int i=0;i(newnumberofvertices); my_elements=xNew(newnumberofelements); @@ -3475,14 +3475,14 @@ /*Start figuring out, out of the partition, which elements belong to this cpu: */ for(int i=0;i(epart); - xDelete(npart); + xDelete(npart); xDelete(index); } /*}}}*/ void FemModel::SmoothedDeviatoricStressTensor(IssmDouble** ptauxx,IssmDouble** ptauyy,IssmDouble** ptauxy){/*{{{*/ - + int elementswidth = this->GetElementsWidth();//just 2D mesh, tria elements int numberofvertices = this->vertices->NumberOfVertices(); IssmDouble weight = 0.; - IssmDouble* tauxx = NULL; - IssmDouble* tauyy = NULL; - IssmDouble* tauxy = NULL; + IssmDouble* tauxx = NULL; + IssmDouble* tauyy = NULL; + IssmDouble* tauxy = NULL; IssmDouble* totalweight = NULL; IssmDouble* deviatoricstressxx = xNew(elementswidth); IssmDouble* deviatoricstressyy = xNew(elementswidth); @@ -3515,10 +3515,10 @@ Vector* vectauyy = new Vector(numberofvertices); Vector* vectauxy = new Vector(numberofvertices); Vector* vectotalweight = new Vector(numberofvertices); - + /*Update the Deviatoric Stress tensor over the elements*/ this->DeviatoricStressx(); - + /*Calculate the Smoothed Deviatoric Stress tensor*/ for(int i=0;ielements->Size();i++){ Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); @@ -3563,7 +3563,7 @@ /*Divide for the total weight*/ for(int i=0;i0); + _assert_(totalweight[i]>0); tauxx[i] = tauxx[i]/totalweight[i]; tauyy[i] = tauyy[i]/totalweight[i]; tauxy[i] = tauxy[i]/totalweight[i]; @@ -3588,7 +3588,7 @@ /*}}}*/ void FemModel::ZZErrorEstimator(IssmDouble** pelementerror){/*{{{*/ - /*Compute the Zienkiewicz and Zhu (ZZ) error estimator for the deviatoric stress tensor. + /*Compute the Zienkiewicz and Zhu (ZZ) error estimator for the deviatoric stress tensor. * Ref.: Zienkiewicz and Zhu, A Simple Error Estimator and Adaptive Procedure for Practical Engineering Analysis, Int. J. Numer. Meth. Eng, 1987*/ IssmDouble Jdet,error,ftxx,ftyy,ftxy; @@ -3608,7 +3608,7 @@ /*Get smoothed deviatoric stress tensor*/ this->SmoothedDeviatoricStressTensor(&smoothedtauxx,&smoothedtauyy,&smoothedtauxy); - + /*Integrate the error over elements*/ for(int i=0;ielements->Size();i++){ Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); @@ -3616,7 +3616,7 @@ element->GetInputListOnVertices(tauyy,DeviatoricStressyyEnum); element->GetInputListOnVertices(tauxy,DeviatoricStressxyEnum); element->GetVerticesSidList(elem_vertices); - + /*Integrate*/ element->GetVerticesCoordinates(&xyz_list); Gauss* gauss=element->NewGauss(2); @@ -3631,9 +3631,9 @@ ftyy+=(tauyy[n]-smoothedtauyy[elem_vertices[n]])*basis[n]; ftxy+=(tauxy[n]-smoothedtauxy[elem_vertices[n]])*basis[n]; } - error+=Jdet*gauss->weight*( pow(ftxx,2)+pow(ftyy,2)+pow(ftxy,2) ); //e^2 + error+=Jdet*gauss->weight*( pow(ftxx,2)+pow(ftyy,2)+pow(ftxy,2) ); //e^2 } - /*Set the error in the global vector*/ + /*Set the error in the global vector*/ sid=element->Sid(); error = sqrt(error);//sqrt(e^2) velementerror->SetValue(sid,error,INS_VAL); @@ -3647,7 +3647,7 @@ /*Serialize and set output*/ (*pelementerror)=velementerror->ToMPISerial(); - + /*Cleanup*/ xDelete(smoothedtauxx); xDelete(smoothedtauyy); @@ -3691,7 +3691,7 @@ /*weight to calculate the smoothed grad H*/ Tria* triaelement = xDynamicCast(element); weight = triaelement->GetArea();//the tria area is a choice for the weight - + /*dH/dx*/ vecdHdx->SetValue(elem_vertices[0],weight*GradH[0],ADD_VAL); vecdHdx->SetValue(elem_vertices[1],weight*GradH[0],ADD_VAL); @@ -3759,7 +3759,7 @@ /*Get smoothed deviatoric stress tensor*/ this->SmoothedGradThickness(&smoothed_dHdx,&smoothed_dHdy); - + /*Integrate the error over elements*/ for(int i=0;ielements->Size();i++){ Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); @@ -3845,7 +3845,7 @@ Vector* vyc = new Vector(numberofelements); IssmDouble* x = NULL; IssmDouble* y = NULL; - IssmDouble* z = NULL; + IssmDouble* z = NULL; IssmDouble* xyz_list = NULL; IssmDouble x1,y1,x2,y2,x3,y3; @@ -3854,7 +3854,7 @@ Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); //element->GetVerticesSidList(elem_vertices); int sid = element->Sid(); - element->GetVerticesCoordinates(&xyz_list); + element->GetVerticesCoordinates(&xyz_list); x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1]; x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1]; x3 = xyz_list[3*2+0];y3 = xyz_list[3*2+1]; @@ -3887,11 +3887,11 @@ if(levelset_type!=MaskGroundediceLevelsetEnum && levelset_type!=MaskIceLevelsetEnum){ _error_("level set type not implemented yet!"); } - + /*Outputs*/ IssmDouble* zerolevelset_points = NULL; int npoints = 0; - + /*Intermediaries*/ int elementswidth = this->GetElementsWidth(); int numberofelements = this->elements->NumberOfElements(); @@ -3904,7 +3904,7 @@ IssmDouble* y_zerolevelset = NULL; int count,sid; IssmDouble xc,yc,x1,y1,x2,y2,x3,y3; - + /*Use the element center coordinate if level set is zero (grounding line or ice front), otherwise set NAN*/ for(int i=0;ielements->Size();i++){ Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); @@ -3911,15 +3911,15 @@ element->GetInputListOnVertices(levelset,levelset_type); element->GetVerticesSidList(elem_vertices); sid= element->Sid(); - element->GetVerticesCoordinates(&xyz_list); + element->GetVerticesCoordinates(&xyz_list); x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1]; x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1]; x3 = xyz_list[3*2+0];y3 = xyz_list[3*2+1]; xc = NAN; - yc = NAN; + yc = NAN; Tria* tria = xDynamicCast(element); if(tria->IsIceInElement()){/*verify if there is ice in the element*/ - if(levelset[0]*levelset[1]<0. || levelset[0]*levelset[2]<0. || + if(levelset[0]*levelset[1]<0. || levelset[0]*levelset[2]<0. || abs(levelset[0]*levelset[1])* pUp, Vector* pNorth, Vector* pEast, Vector* pX, Vector* pY, IssmDouble* xx, IssmDouble* yy){/*{{{*/ int ns,nsmax; - + /*Go through elements, and add contribution from each element to the deflection vector wg:*/ ns = elements->Size(); - + /*Figure out max of ns: */ ISSM_MPI_Reduce(&ns,&nsmax,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm()); ISSM_MPI_Bcast(&nsmax,1,ISSM_MPI_INT,0,IssmComm::GetComm()); @@ -4104,7 +4104,7 @@ pY->Assemble(); } } - + /*One last time: */ pUp->Assemble(); pNorth->Assemble(); @@ -4123,11 +4123,11 @@ IssmDouble eartharea_cpu=0; int ns,nsmax; - + /*Go through elements, and add contribution from each element to the deflection vector wg:*/ ns = elements->Size(); - - /*First, figure out the surface area of Earth: */ + + /*First, figure out the surface area of Earth: */ for(int i=0;i(elements->GetObjectByOffset(i)); eartharea_cpu += element->GetAreaSpherical(); @@ -4151,7 +4151,7 @@ pEast->Assemble(); } } - + /*One last time: */ pUp->Assemble(); pNorth->Assemble(); @@ -4213,7 +4213,7 @@ if(i%loop==0)pRSLgi->Assemble(); } if(VerboseConvergence())_printf0_("\n"); - + /*One last time: */ pRSLgi->Assemble(); @@ -4231,12 +4231,12 @@ /*serialized vectors:*/ IssmDouble* RSLg_old=NULL; - + IssmDouble eartharea=0; IssmDouble eartharea_cpu=0; int ns,nsmax; - + /*Serialize vectors from previous iteration:*/ RSLg_old=pRSLg_old->ToMPISerial(); @@ -4248,7 +4248,7 @@ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); eartharea_cpu += element->GetAreaSpherical(); } - + ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); @@ -4266,7 +4266,7 @@ if(i%loop==0)pRSLgo->Assemble(); } if(verboseconvolution)if(VerboseConvergence())_printf0_("\n"); - + /*Free ressources:*/ xDelete(RSLg_old); } @@ -4354,8 +4354,8 @@ /*Free ressources:*/ xDelete(RSLg_old); - - + + } /*}}}*/ void FemModel::SealevelriseElastic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pRSLg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/ @@ -4362,18 +4362,18 @@ /*serialized vectors:*/ IssmDouble* RSLg=NULL; - + IssmDouble eartharea=0; IssmDouble eartharea_cpu=0; int ns,nsmax; - + /*Serialize vectors from previous iteration:*/ RSLg=pRSLg->ToMPISerial(); /*Go through elements, and add contribution from each element to the deflection vector wg:*/ ns = elements->Size(); - + /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */ for(int i=0;i(elements->GetObjectByOffset(i)); @@ -4401,7 +4401,7 @@ } } } - + /*One last time: */ pUp->Assemble(); if(horiz){ @@ -4435,13 +4435,13 @@ } ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - + ISSM_MPI_Reduce (&oceanvalue_cpu,&oceanvalue,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); ISSM_MPI_Bcast(&oceanvalue,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); /*Free ressources:*/ xDelete(RSLg_serial); - + return oceanvalue/oceanarea; } /*}}}*/ @@ -4457,7 +4457,7 @@ IssmDouble* old_active = NULL; int* eplzigzag_counter = NULL; int eplflip_lock; - + HydrologyDCEfficientAnalysis* effanalysis = new HydrologyDCEfficientAnalysis(); HydrologyDCInefficientAnalysis* inefanalysis = new HydrologyDCInefficientAnalysis(); @@ -4464,10 +4464,10 @@ /*Step 1: update mask, the mask might be extended by residual and/or using downstream sediment head*/ mask=new Vector(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum)); recurence=new Vector(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum)); - this->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum); - this->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum); + this->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum); + this->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum); GetVectorFromInputsx(&old_active,this,HydrologydcMaskEplactiveNodeEnum,NodeSIdEnum); - + for (int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); effanalysis->HydrologyEPLGetMask(mask,recurence,eplzigzag_counter,element); @@ -4485,8 +4485,8 @@ this->parameters->SetParam(eplzigzag_counter,this->nodes->Size(),EplZigZagCounterEnum); /*Assemble and serialize*/ mask->Assemble(); - serial_mask=mask->ToMPISerial(); - + serial_mask=mask->ToMPISerial(); + xDelete(eplzigzag_counter); xDelete(serial_rec); xDelete(old_active); @@ -4528,7 +4528,7 @@ delete inefanalysis; int sum_counter; ISSM_MPI_Reduce(&counter,&sum_counter,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm()); counter=sum_counter; *pEplcount = counter; if(VerboseSolution()) _printf0_(" Number of active nodes in EPL layer: "<< counter <<"\n"); @@ -4649,7 +4649,7 @@ xDelete(serial_active); int sum_counter; ISSM_MPI_Reduce(&counter,&sum_counter,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm()); counter=sum_counter; *pL2count = counter; if(VerboseSolution()) _printf0_(" Number of active nodes L2 Projection: "<< counter <<"\n"); @@ -4734,7 +4734,7 @@ } } /*}}}*/ -#ifdef _HAVE_JAVASCRIPT_ +#ifdef _HAVE_JAVASCRIPT_ FemModel::FemModel(IssmDouble* buffer, int buffersize, char* toolkits, char* solution, char* modelname,ISSM_MPI_Comm incomm, bool trace){ /*{{{*/ /*configuration: */ int solution_type; @@ -4749,12 +4749,12 @@ /*From command line arguments, retrieve different filenames needed to create the FemModel: */ solution_type=StringToEnumx(solution); - + /*Create femmodel from input files: */ profiler->Start(MPROCESSOR); this->InitFromBuffers((char*)buffer,buffersize,toolkits, solution_type,trace,NULL); profiler->Stop(MPROCESSOR); - + /*Save communicator in the parameters dataset: */ this->parameters->AddObject(new GenericParam(incomm,FemModelCommEnum)); @@ -4769,7 +4769,7 @@ char** poutputbuffer; size_t* poutputbuffersize; - + /*Before we delete the profiler, report statistics for this run: */ profiler->Stop(TOTAL); //final tagging _printf0_("\n"); @@ -4782,7 +4782,7 @@ <TotalTimeModSec(TOTAL)<<" sec" ); _printf0_("\n"); - + /*Before we close the output file, recover the buffer and size:*/ outputbufferparam = xDynamicCast*>(this->parameters->FindParamObject(OutputBufferPointerEnum)); poutputbuffer=outputbufferparam->GetParameterValue(); @@ -4822,8 +4822,8 @@ fclose(toolkitsoptionsfid); /*Open output file once for all and add output file descriptor to parameters*/ - output_fid=open_memstream(&outputbuffer,&outputsize); - if(output_fid==NULL||output_fid==0)_error_("could not initialize output stream"); + output_fid=open_memstream(&outputbuffer,&outputsize); + if(output_fid==NULL)_error_("could not initialize output stream"); this->parameters->SetParam(output_fid,OutputFilePointerEnum); this->parameters->AddObject(new GenericParam(&outputbuffer,OutputBufferPointerEnum)); this->parameters->AddObject(new GenericParam(&outputsize,OutputBufferSizePointerEnum)); @@ -4864,7 +4864,7 @@ this->amrbamg->thicknesserror_threshold>0||this->amrbamg->deviatoricerror_threshold>0){ /*Initialize hmaxvertices with NAN*/ hmaxvertices_serial=xNew(numberofvertices); - for(int i=0;iamrbamg->groundingline_distance>0) this->GethmaxVerticesFromZeroLevelSetDistance(hmaxvertices_serial,MaskGroundediceLevelsetEnum); if(this->amrbamg->icefront_distance>0) this->GethmaxVerticesFromZeroLevelSetDistance(hmaxvertices_serial,MaskIceLevelsetEnum); @@ -4871,7 +4871,7 @@ if(this->amrbamg->thicknesserror_threshold>0) this->GethmaxVerticesFromEstimators(hmaxvertices_serial,ThicknessErrorEstimatorEnum); if(this->amrbamg->deviatoricerror_threshold>0) this->GethmaxVerticesFromEstimators(hmaxvertices_serial,DeviatoricStressErrorEstimatorEnum); } - + if(my_rank==0){ this->amrbamg->ExecuteRefinementBamg(vector_serial,hmaxvertices_serial,&newnumberofvertices,&newnumberofelements,&newx,&newy,&newz,&newelementslist); if(newnumberofvertices<=0 || newnumberofelements<=0) _error_("Error in the refinement process."); @@ -4891,10 +4891,10 @@ newz=xNew(newnumberofvertices); newelementslist=xNew(newnumberofelements*this->GetElementsWidth()); } - ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm()); /*Assign output pointers*/ *pnewnumberofvertices = newnumberofvertices; @@ -4927,7 +4927,7 @@ /*Create bamg data structures for bamg*/ this->amrbamg = new AmrBamg(); - + /*Get amr parameters*/ this->parameters->FindParam(&hmin,AmrHminEnum); this->parameters->FindParam(&hmax,AmrHmaxEnum); @@ -4951,7 +4951,7 @@ this->amrbamg->SetBamgOpts(hmin,hmax,err,gradation); /*Re-create original mesh and put it in bamg structure (only cpu 0)*/ - if(my_rank==0){ + if(my_rank==0){ this->amrbamg->Initialize(elements,x,y,numberofvertices,numberofelements); } @@ -4965,7 +4965,7 @@ void FemModel::GethmaxVerticesFromZeroLevelSetDistance(IssmDouble* hmaxvertices,int levelset_type){/*{{{*/ if(!hmaxvertices) _error_("hmaxvertices is NULL!\n"); - + /*Intermediaries*/ int numberofvertices = this->vertices->NumberOfVertices(); IssmDouble* verticedistance = NULL; @@ -4972,7 +4972,7 @@ IssmDouble threshold,resolution; switch(levelset_type){ - case MaskGroundediceLevelsetEnum: + case MaskGroundediceLevelsetEnum: threshold = this->amrbamg->groundingline_distance; resolution = this->amrbamg->groundingline_resolution; break; @@ -4986,7 +4986,7 @@ /*Get vertice distance to zero level set points*/ this->GetVerticeDistanceToZeroLevelSet(&verticedistance,levelset_type); if(!verticedistance) _error_("verticedistance is NULL!\n"); - + /*Fill hmaxVertices*/ for(int i=0;ielements->NumberOfElements(); int numberofvertices = this->vertices->NumberOfVertices(); IssmDouble* maxlength = xNew(numberofelements); - IssmDouble* error_vertices = xNewZeroInit(numberofvertices); + IssmDouble* error_vertices = xNewZeroInit(numberofvertices); IssmDouble* error_elements = NULL; IssmDouble* x = NULL; IssmDouble* y = NULL; @@ -5021,7 +5021,7 @@ /*Fill variables*/ switch(errorestimator_type){ - case ThicknessErrorEstimatorEnum: + case ThicknessErrorEstimatorEnum: threshold = this->amrbamg->thicknesserror_threshold; groupthreshold = this->amrbamg->thicknesserror_groupthreshold; resolution = this->amrbamg->thicknesserror_resolution; @@ -5046,12 +5046,12 @@ switch(errorestimator_type){ case ThicknessErrorEstimatorEnum: this->amrbamg->thicknesserror_maximum = maxerror;break; case DeviatoricStressErrorEstimatorEnum: this->amrbamg->deviatoricerror_maximum = maxerror;break; - } + } } /*Get mesh*/ this->GetMesh(this->vertices,this->elements,&x,&y,&z,&index); - + /*Fill error_vertices (this is the sum of all elements connected to the vertex)*/ for(int i=0;igroupthreshold*maxerror) refine=true; } } - /*Now, fill the hmaxvertices if requested*/ + /*Now, fill the hmaxvertices if requested*/ if(refine){ for(int j=0;jvertices->NumberOfVertices(); IssmDouble* levelset_points= NULL; @@ -5121,8 +5121,8 @@ /*Get vertices coordinates*/ VertexCoordinatesx(&x,&y,&z,this->vertices,false) ; - - /*Get points which level set is zero (center of elements with zero level set)*/ + + /*Get points which level set is zero (center of elements with zero level set)*/ this->GetZeroLevelSetPoints(&levelset_points,numberofpoints,levelset_type); /*Find the minimal vertice distance to the zero levelset (grounding line or ice front)*/ @@ -5131,9 +5131,9 @@ verticedistance[i]=INFINITY; for(int j=0;jamr->groundingline_distance>0) this->GetElementDistanceToZeroLevelSet(&gl_distance,MaskGroundediceLevelsetEnum); if(this->amr->icefront_distance>0) this->GetElementDistanceToZeroLevelSet(&if_distance,MaskIceLevelsetEnum); - if(this->amr->thicknesserror_threshold>0) this->ThicknessZZErrorEstimator(&thicknesserror); - if(this->amr->deviatoricerror_threshold>0) this->ZZErrorEstimator(&deviatoricerror); - + if(this->amr->thicknesserror_threshold>0) this->ThicknessZZErrorEstimator(&thicknesserror); + if(this->amr->deviatoricerror_threshold>0) this->ZZErrorEstimator(&deviatoricerror); + if(my_rank==0){ this->amr->ExecuteRefinement(gl_distance,if_distance,deviatoricerror,thicknesserror, - &newnumberofvertices,&newnumberofelements,&newx,&newy,&newelementslist); + &newnumberofvertices,&newnumberofelements,&newx,&newy,&newelementslist); newz=xNewZeroInit(newnumberofvertices); if(newnumberofvertices<=0 || newnumberofelements<=0) _error_("Error in the ReMeshNeopz."); } @@ -5186,12 +5186,12 @@ newz=xNew(newnumberofvertices); newelementslist=xNew(newnumberofelements*this->GetElementsWidth()); } - ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm()); - /*Assign the pointers*/ + /*Assign the pointers*/ (*pnewelementslist) = newelementslist; //Matlab indexing (*pnewx) = newx; (*pnewy) = newy; @@ -5207,7 +5207,7 @@ } /*}}}*/ void FemModel::InitializeAdaptiveRefinementNeopz(void){/*{{{*/ - + /*Define variables*/ int my_rank = IssmComm::GetRank(); int numberofvertices = this->vertices->NumberOfVertices(); @@ -5216,7 +5216,7 @@ IssmDouble* y = NULL; IssmDouble* z = NULL; int* elements = NULL; - + /*Initialize field as NULL for now*/ this->amr = NULL; @@ -5224,12 +5224,12 @@ /*elements comes in Matlab indexing*/ this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elements); - /*Create initial mesh (coarse mesh) in neopz data structure*/ + /*Create initial mesh (coarse mesh) in neopz data structure*/ /*Just CPU #0 should keep AMR object*/ /*Initialize refinement pattern*/ this->SetRefPatterns(); this->amr = new AdaptiveMeshRefinement(); - this->amr->refinement_type=1;//1 is refpattern; 0 is uniform (faster) + this->amr->refinement_type=1;//1 is refpattern; 0 is uniform (faster) /*Get amr parameters*/ this->parameters->FindParam(&this->amr->level_max,AmrLevelMaxEnum); this->parameters->FindParam(&this->amr->gradation,AmrGradationEnum); @@ -5242,7 +5242,7 @@ this->parameters->FindParam(&this->amr->deviatoricerror_threshold,AmrDeviatoricErrorThresholdEnum); this->parameters->FindParam(&this->amr->deviatoricerror_groupthreshold,AmrDeviatoricErrorGroupThresholdEnum); this->parameters->FindParam(&this->amr->deviatoricerror_maximum,AmrDeviatoricErrorMaximumEnum); - if(my_rank==0){ + if(my_rank==0){ this->amr->CreateInitialMesh(numberofvertices,numberofelements,x,y,elements); } @@ -5263,7 +5263,7 @@ /*Output*/ IssmDouble* elementdistance; - + /*Intermediaries*/ int numberofelements = this->elements->NumberOfElements(); Vector* velementdistance = new Vector(numberofelements); @@ -5272,14 +5272,14 @@ IssmDouble mindistance,distance; IssmDouble xc,yc,x1,y1,x2,y2,x3,y3; int numberofpoints; - - /*Get points which level set is zero (center of elements with zero level set, levelset_points is serial)*/ + + /*Get points which level set is zero (center of elements with zero level set, levelset_points is serial)*/ this->GetZeroLevelSetPoints(&levelset_points,numberofpoints,levelset_type); for(int i=0;ielements->Size();i++){//parallel Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); int sid = element->Sid(); - element->GetVerticesCoordinates(&xyz_list); + element->GetVerticesCoordinates(&xyz_list); x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1]; x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1]; x3 = xyz_list[3*2+0];y3 = xyz_list[3*2+1]; @@ -5289,11 +5289,11 @@ /*Loop over each point (where level set is zero)*/ for(int j=0;jSetValue(sid,mindistance,INS_VAL); xDelete(xyz_list); - } + } /*Assemble*/ velementdistance->Assemble();