Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 21801) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 21802) @@ -47,4 +47,5 @@ ./bamg/Mesh.cpp\ ./shared/Bamg/BigPrimeNumber.cpp\ + ./classes/AmrBamg.cpp\ ./modules/Bamgx/Bamgx.cpp\ ./modules/BamgConvertMeshx/BamgConvertMeshx.cpp\ Index: /issm/trunk-jpl/src/c/bamg/BamgOpts.cpp =================================================================== --- /issm/trunk-jpl/src/c/bamg/BamgOpts.cpp (revision 21801) +++ /issm/trunk-jpl/src/c/bamg/BamgOpts.cpp (revision 21802) @@ -5,28 +5,28 @@ BamgOpts::BamgOpts(){/*{{{*/ - this->anisomax=0; - this->cutoff=0; - this->coeff=0; - this->errg=0; - this->gradation=0; - this->Hessiantype=0; - this->MaxCornerAngle=0; - this->maxnbv=0; - this->maxsubdiv=0; - this->Metrictype=0; - this->nbjacobi=0; - this->nbsmooth=0; - this->omega=0; - this->power=0; - this->random=0; - this->verbose=0; + this->anisomax = 0; + this->cutoff = 0; + this->coeff = 0; + this->errg = 0; + this->gradation = 0; + this->Hessiantype = 0; + this->MaxCornerAngle = 0; + this->maxnbv = 0; + this->maxsubdiv = 0; + this->Metrictype = 0; + this->nbjacobi = 0; + this->nbsmooth = 0; + this->omega = 0; + this->power = 0; + this->random = 0; + this->verbose = 0; - this->Crack=0; - this->geometricalmetric=0; - this->KeepVertices=0; - this->splitcorners=0; + this->Crack = 0; + this->geometricalmetric = 0; + this->KeepVertices = 0; + this->splitcorners = 0; - this->hmin=0; - this->hmax=0; + this->hmin = 0; + this->hmax = 0; this->hminVertices=NULL; this->hminVerticesSize[0]=this->hminVerticesSize[1]=0; this->hmaxVertices=NULL; this->hmaxVerticesSize[0]=this->hmaxVerticesSize[1]=0; Index: /issm/trunk-jpl/src/c/bamg/Mesh.cpp =================================================================== --- /issm/trunk-jpl/src/c/bamg/Mesh.cpp (revision 21801) +++ /issm/trunk-jpl/src/c/bamg/Mesh.cpp (revision 21802) @@ -698,5 +698,5 @@ for (i=0;iSubDomains[i*4+0]=3; - bamgmesh->SubDomains[i*4+1]=reft[GetId(subdomains[i].head)]; + bamgmesh->SubDomains[i*4+1]=reft[GetId(subdomains[i].head)]+1;//MATLAB indexing bamgmesh->SubDomains[i*4+2]=1; bamgmesh->SubDomains[i*4+3]=subdomains[i].ReferenceNumber; @@ -3076,8 +3076,8 @@ IssmDouble area_3 =((bv.r.x -(*tcvj)(1)->r.x)*((*tcvj)(0)->r.y-(*tcvj)(1)->r.y) - (bv.r.y -(*tcvj)(1)->r.y)*((*tcvj)(0)->r.x-(*tcvj)(1)->r.x)); - if(area_1<0 || area_2<0 || area_3<0){ - pointsoutside = true; - continue; - } + // if(area_1<0 || area_2<0 || area_3<0){ + // pointsoutside = true; + // continue; + // } if(!bv.GeomEdgeHook){ vertices[nbv].r = bv.r; @@ -3086,5 +3086,5 @@ } } - if(pointsoutside) _printf_("WARNING: One or more points of the initial mesh fall outside of the geometric boundary\n"); + //if(pointsoutside) _printf_("WARNING: One or more points of the initial mesh fall outside of the geometric boundary\n"); Bh.CreateSingleVertexToTriangleConnectivity(); InsertNewPoints(nbvold,NbTSwap,bamgopts->random); Index: /issm/trunk-jpl/src/c/classes/AmrBamg.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/AmrBamg.cpp (revision 21802) +++ /issm/trunk-jpl/src/c/classes/AmrBamg.cpp (revision 21802) @@ -0,0 +1,131 @@ +/*!\file AmrBamg.cpp + * \brief: implementation of the adaptive mesh refinement tool based on bamg + */ + +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "./AmrBamg.h" +#include "../bamg/bamgobjects.h" +#include "../modules/Bamgx/Bamgx.h" + +using namespace bamg; +using namespace std; + +/*Constructor, copy, clean up and destructor*/ +AmrBamg::AmrBamg(IssmDouble hmin, IssmDouble hmax){/*{{{*/ + + this->geometry = NULL; + this->fathermesh = NULL; + this->previousmesh = NULL; + + /*Only initialize options for now (same as bamg.m)*/ + this->options = new BamgOpts(); + this->options->anisomax = 10.e30; + this->options->cutoff = 10.e-5; + this->options->coeff = 1; + this->options->errg = 0.1; + this->options->gradation = 1.5; + this->options->Hessiantype = 0; + this->options->MaxCornerAngle = 10; + this->options->maxnbv = 1e6; + this->options->maxsubdiv = 10; + this->options->Metrictype = 0; + this->options->nbjacobi = 1; + this->options->nbsmooth = 3; + this->options->omega = 1.8; + this->options->power = 1; + this->options->random = 0; + this->options->verbose = 0; + this->options->Crack = 0; + this->options->geometricalmetric = 0; + this->options->KeepVertices = 1; /*!!!!! VERY IMPORTANT !!!!!*/ + this->options->splitcorners = 1; + this->options->hmin = hmin; + this->options->hmax = hmax; +} +/*}}}*/ +AmrBamg::~AmrBamg(){/*{{{*/ + + if(this->geometry) delete this->geometry; + if(this->fathermesh) delete this->fathermesh; + if(this->previousmesh) delete this->previousmesh; + if(this->options) delete this->options; + +} +/*}}}*/ + +/*Methods*/ +void AmrBamg::Initialize(int* elements,IssmDouble* x,IssmDouble* y,int numberofvertices,int numberofelements){/*{{{*/ + + /*Check options*/ + _assert_(this->options); + this->options->Check(); + + /*Read father mesh and create geometry*/ + Mesh* Th=new Mesh(elements,x,y,numberofvertices,numberofelements); + + /*Write geometry*/ + this->geometry = new BamgGeom(); + Th->Gh.WriteGeometry(this->geometry,this->options); + + /*Write father mesh*/ + this->fathermesh = new BamgMesh(); + Th->WriteMesh(this->fathermesh,this->options); + + /*Cleanup and return*/ + delete Th; +}/*}}}*/ +void AmrBamg::ExecuteRefinementBamg(int* pnewnumberofvertices,int *pnewnumberofelements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist){/*{{{*/ + + /*Intermediaries*/ + BamgGeom* geomout=new BamgGeom(); + BamgMesh* meshout=new BamgMesh(); + + /*Some checks*/ + _assert_(this->geometry); + _assert_(this->options); + _assert_(this->fathermesh); + + /*remesh*/ + if(this->previousmesh){ + Bamgx(meshout,geomout,this->previousmesh,this->geometry,this->options); + } + else{ + Bamgx(meshout,geomout,this->fathermesh,this->geometry,this->options); + } + + /*Change previous mesh*/ + if(this->previousmesh) delete this->previousmesh; + this->previousmesh = meshout; + + /*Prepare output*/ + int nbv = meshout->VerticesSize[0]; + int nbt = meshout->TrianglesSize[0]; + IssmDouble *x = xNew(nbv); + IssmDouble *y = xNew(nbv); + IssmDouble *z = xNew(nbv); + for(int i=0;iVertices[i*2+0]; + y[i] = meshout->Vertices[i*2+1]; + z[i] = 0.; + } + int* elementslist= xNew(nbt*3); + for(int i=0;i(meshout->Triangles[4*i+0]); + elementslist[3*i+1] = reCast(meshout->Triangles[4*i+1]); + elementslist[3*i+2] = reCast(meshout->Triangles[4*i+2]); + } + + /*Cleanup and return*/ + delete geomout; + *pnewnumberofvertices = nbv; + *pnewnumberofelements = nbt; + *px = x; + *py = y; + *pz = z; + *pelementslist = elementslist; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/AmrBamg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/AmrBamg.h (revision 21802) +++ /issm/trunk-jpl/src/c/classes/AmrBamg.h (revision 21802) @@ -0,0 +1,30 @@ +#ifndef AMRBAMG +#define AMRBAMG + +/*Includes*/ +#include "../shared/shared.h" +#include "../toolkits/toolkits.h" + +#include "../bamg/BamgMesh.h" +#include "../bamg/BamgGeom.h" +#include "../bamg/BamgOpts.h" + +class AmrBamg{ + + public: + /* Constructor, destructor etc*/ + AmrBamg(IssmDouble hmin, IssmDouble hmax); + ~AmrBamg(); + + /*General methods*/ + void Initialize(int* elements,IssmDouble* x,IssmDouble* y,int numberofvertices,int numberofelements); + void ExecuteRefinementBamg(int* pnewnumberofvertices,int *pnewnumberofelements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist); + + private: + BamgGeom* geometry; + BamgMesh* fathermesh; + BamgMesh* previousmesh; + BamgOpts* options; +}; + +#endif Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 21801) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 21802) @@ -45,5 +45,5 @@ /*configuration: */ - int solution_type; + int solution_type,amrtype,amr_frequency; int ierr; @@ -80,7 +80,21 @@ this->parameters->AddObject(new GenericParam(incomm,FemModelCommEnum)); - #ifdef _HAVE_NEOPZ_ - this->InitializeAdaptiveRefinement(); - #endif + /*AMR stuff*/ + this->parameters->FindParam(&amr_frequency,TransientAmrFrequencyEnum); + if(amr_frequency){ + this->parameters->FindParam(&amrtype,AmrTypeEnum); + switch(amrtype){ + + #ifdef _HAVE_NEOPZ_ + case AmrNeopzEnum: this->InitializeAdaptiveRefinementNeopz(); break; + #endif + + #ifdef _HAVE_BAMG_ + case AmrBamgEnum: this->InitializeAdaptiveRefinementBamg(); break; + #endif + + default: _error_("not implemented yet"); + } + } /*Free resources */ @@ -140,4 +154,8 @@ #ifdef _HAVE_NEOPZ_ if(amr)delete amr; + #endif + + #ifdef _HAVE_BAMG_ + if(amrbamg)delete amrbamg; #endif @@ -1502,4 +1520,148 @@ } /*}}}*/ +void FemModel::ReMesh(void){/*{{{*/ + + /*Intermediaries*/ + IssmDouble *newx = NULL; + IssmDouble *newy = NULL; + IssmDouble *newz = NULL; + int *newelementslist = NULL; + int newnumberofvertices = -1; + int newnumberofelements = -1; + bool* my_elements = NULL; + int* my_vertices = NULL; + int elementswidth = this->GetElementsWidth();//just tria elements in this version + int amrtype; + + /*Branch to specific amr depending on requested method*/ + parameters->FindParam(&amrtype,AmrTypeEnum); + switch(amrtype){ + #ifdef _HAVE_NEOPZ_ + case AmrNeopzEnum: this->ReMeshNeopz(&newnumberofvertices,&newnumberofelements,&newx,&newy,&newz,&newelementslist); break; + #endif + + #ifdef _HAVE_BAMG_ + case AmrBamgEnum: this->ReMeshBamg(&newnumberofvertices,&newnumberofelements,&newx,&newy,&newz,&newelementslist); break; + #endif + + 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); + + if(this->loads->Size()!=0) _error_("not supported yet"); + + /*Create vertices*/ + Vertices* new_vertices=new Vertices(); + this->CreateVertices(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,my_vertices,newx,newy,newz,new_vertices); + + /*Creating elements*/ + /*Just Tria in this version*/ + Elements* new_elements=new Elements(); + this->CreateElements(newnumberofelements,elementswidth,newelementslist,my_elements,new_elements); + + /*Creating materials*/ + Materials* new_materials=new Materials(); + this->CreateMaterials(newnumberofelements,my_elements,new_materials); + + /*Creating nodes and constraints*/ + /*Just SSA (2D) and P1 in this version*/ + Nodes* new_nodes=new Nodes(); + Constraints* new_constraints=new Constraints(); + + int nodecounter=0; + int constraintcounter=0; + for(int i=0;inummodels;i++){//create nodes for each analysis in analysis_type_list + + int analysis_enum = this->analysis_type_list[i]; + + /*As the domain is 2D, it is not necessary to create nodes for this analysis*/ + /*itapopo must verify if domain is not 3D. Only 2D in this version!*/ + if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue; + + this->CreateNodes(newnumberofvertices,my_vertices,nodecounter,analysis_enum,new_nodes); + if(analysis_enum==StressbalanceAnalysisEnum) this->CreateConstraints(newnumberofvertices,newnumberofelements,nodecounter,constraintcounter,newx,newy,my_vertices,new_constraints); + this->UpdateElements(newnumberofelements,newelementslist,my_elements,nodecounter,i,new_elements); + + if(new_nodes->Size()) nodecounter=new_nodes->MaximumId(); + constraintcounter = new_constraints->NumberOfConstraints(); + /*Make sure nodecounter is at least 0 (if no node exists, maxid will be -1*/ + _assert_(nodecounter>=0); + } + + new_elements->Presort(); + new_nodes->Presort(); + new_vertices->Presort(); + this->loads->Presort(); + new_materials->Presort(); + new_constraints->Presort(); + + /*reset hooks for elements, loads and nodes: */ + new_elements->ResetHooks(); + this->loads->ResetHooks(); + new_materials->ResetHooks(); + + /*do the post-processing of the datasets to get an FemModel that can actually run analyses: */ + int analysis_type; + for(int i=0;inummodels;i++){ + analysis_type=this->analysis_type_list[i]; + //SetCurrentConfiguration(analysis_type); + + 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); + this->parameters->SetParam(analysis_type,ConfigurationTypeEnum); + + /*configure elements, loads and nodes, for this new analysis: */ + new_elements->SetCurrentConfiguration(new_elements,this->loads,new_nodes,new_vertices,new_materials,this->parameters); + this->loads->SetCurrentConfiguration(new_elements,this->loads,new_nodes,new_vertices,new_materials,this->parameters); + + /*take care of toolkits options, that depend on this analysis type (present only after model processor)*/ + if(this->parameters->Exist(ToolkitsOptionsStringsEnum)){ + ToolkitsOptionsFromAnalysis(this->parameters,analysis_type); + if(VerboseSolver()) _printf0_(" toolkits Options set for analysis type: " << EnumToStringx(analysis_type) << "\n"); + } + + ConfigureObjectsx(new_elements,this->loads,new_nodes,new_vertices,new_materials,this->parameters); + 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); + NodesDofx(new_nodes,this->parameters,analysis_type); + } + + /*Finally: interpolate all inputs and insert them into the new elements.*/ + this->InterpolateInputs(new_vertices,new_elements); + + /*Delete old structure and set new pointers*/ + delete this->vertices; this->vertices = new_vertices; + delete this->elements; this->elements = new_elements; + delete this->nodes; this->nodes = new_nodes; + delete this->constraints; this->constraints = new_constraints; + delete this->materials; this->materials = new_materials; + + GetMaskOfIceVerticesLSMx(this); + + /*Insert MISMIP+ bed topography*/ + if(false) this->BedrockFromMismipPlus(); + + /*Adjust base, thickness and mask grounded ice leve set*/ + if(true) this->AdjustBaseThicknessAndMask(); + + /*Reset current configuration: */ + analysis_type=this->analysis_type_list[this->analysis_counter]; + SetCurrentConfiguration(analysis_type); + + /*Cleanup*/ + xDelete(newx); + xDelete(newy); + xDelete(newz); + xDelete(newelementslist); + xDelete(my_vertices); + xDelete(my_elements); +} +/*}}}*/ void FemModel::RequestedDependentsx(void){/*{{{*/ @@ -2234,947 +2396,6 @@ } /*}}}*/ -#ifdef _HAVE_DAKOTA_ -void FemModel::DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses){/*{{{*/ - - int i,j; - int my_rank; - - /*intermediary: */ - char root[50]; - int index; - int npart; - double femmodel_response; - int flag; - double *vertex_response = NULL; - double *qmu_response = NULL; - double *responses_pointer = NULL; - - /*retrieve npart: */ - parameters->FindParam(&npart,QmuNumberofpartitionsEnum); - my_rank=IssmComm::GetRank(); - - /*save the d_responses pointer: */ - responses_pointer=d_responses; - - //watch out, we have more d_numresponses than numresponsedescriptors, because the responses have been expanded if they were scaled. - //because we don't know the d_responses descriptors (the scaled ones) we can't key off them, so we will key off the responses_descriptors: */ - - for(i=0;i(vertex_response); - xDelete(qmu_response); - - } - else if (flag==IndexedEnum){ - - /*indexed response: plug index into parameters and call response module: */ - parameters->SetParam(index,IndexEnum); - - this->Responsex(&femmodel_response,root); - - if(my_rank==0){ - /*plug response: */ - responses_pointer[0]=femmodel_response; - - /*increment response_pointer :*/ - responses_pointer++; - } - } - else if (flag==NodalEnum){ - _error_("nodal response functions not supported yet!"); - - /*increment response_pointer :*/ - responses_pointer++; - } - else if (flag==RegularEnum){ - - /*perfectly normal response function: */ - this->Responsex(&femmodel_response,root); - - if(my_rank==0){ - /*plug response: */ - responses_pointer[0]=femmodel_response; - - /*increment response_pointer :*/ - responses_pointer++; - } - } - else _error_("flag type " << flag << " not supported yet for response analysis"); - } - - /*Synthesize echo: {{{*/ - if(my_rank==0){ - _printf_(" responses: " << d_numresponses << ": "); - for(i=0;i* wg,Vector* dwgdt, IssmDouble* x, IssmDouble* y){ /*{{{*/ - - /*Go through elements, and add contribution from each element to the deflection vector wg:*/ - for(int i=0;iSize();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->GiaDeflection(wg,dwgdt, x,y); - } -} -/*}}}*/ -#endif -#ifdef _HAVE_ESA_ -void FemModel::EsaGeodetic2D(Vector* pUp, Vector* pNorth, Vector* pEast, 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()); - - /*Call the esa geodetic core: */ - for(int i=0;i(elements->GetObjectByOffset(i)); - element->EsaGeodetic2D(pUp,pNorth,pEast,xx,yy); - } - if(i%100==0){ - pUp->Assemble(); - pNorth->Assemble(); - pEast->Assemble(); - } - } - - /*One last time: */ - pUp->Assemble(); - pNorth->Assemble(); - pEast->Assemble(); - - /*Free ressources:*/ - xDelete(xx); - xDelete(yy); -} -/*}}}*/ -void FemModel::EsaGeodetic3D(Vector* pUp, Vector* pNorth, Vector* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){/*{{{*/ - - IssmDouble eartharea=0; - 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: */ - for(int i=0;i(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()); - - /*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()); - - /*Call the esa geodetic core: */ - for(int i=0;i(elements->GetObjectByOffset(i)); - element->EsaGeodetic3D(pUp,pNorth,pEast,latitude,longitude,radius,xx,yy,zz,eartharea); - } - if(i%100==0){ - pUp->Assemble(); - pNorth->Assemble(); - pEast->Assemble(); - } - } - - /*One last time: */ - pUp->Assemble(); - pNorth->Assemble(); - pEast->Assemble(); - - /*Free ressources:*/ - xDelete(latitude); - xDelete(longitude); - xDelete(radius); - xDelete(xx); - xDelete(yy); - xDelete(zz); -} -/*}}}*/ -#endif -#ifdef _HAVE_SEALEVELRISE_ -void FemModel::SealevelriseEustatic(Vector* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius) { /*{{{*/ - - /*serialized vectors:*/ - IssmDouble eustatic = 0.; - IssmDouble eustatic_cpu = 0.; - IssmDouble eustatic_cpu_e = 0.; - IssmDouble oceanarea = 0.; - IssmDouble oceanarea_cpu = 0.; - IssmDouble eartharea = 0.; - 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 area of the ocean, which is needed to compute the eustatic component: */ - for(int i=0;i(elements->GetObjectByOffset(i)); - oceanarea_cpu += element->OceanArea(); - eartharea_cpu += element->GetAreaSpherical(); - } - 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()); - _assert_(oceanarea>0.); - - 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()); - - /*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()); - - /*Call the sea level rise core: */ - for(int i=0;i(elements->GetObjectByOffset(i)); - element->SealevelriseEustatic(pSgi,&eustatic_cpu_e,latitude,longitude,radius,oceanarea,eartharea); - eustatic_cpu+=eustatic_cpu_e; - } - if(i%100==0)pSgi->Assemble(); - } - if(VerboseConvergence())_printf0_("\n"); - - /*One last time: */ - pSgi->Assemble(); - - /*Sum all eustatic components from all cpus:*/ - ISSM_MPI_Reduce (&eustatic_cpu,&eustatic,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&eustatic,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - _assert_(!xIsNan(eustatic)); - - /*Assign output pointers:*/ - *peustatic=eustatic; - -} -/*}}}*/ -void FemModel::SealevelriseNonEustatic(Vector* pSgo, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution){/*{{{*/ - - /*serialized vectors:*/ - IssmDouble* Sg_old=NULL; - - IssmDouble eartharea=0; - IssmDouble eartharea_cpu=0; - - int ns,nsmax; - - /*Serialize vectors from previous iteration:*/ - Sg_old=pSg_old->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)); - 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()); - - /*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()); - - /*Call the sea level rise core: */ - for(int i=0;i(elements->GetObjectByOffset(i)); - element->SealevelriseNonEustatic(pSgo,Sg_old,latitude,longitude,radius,eartharea); - } - if(i%100==0)pSgo->Assemble(); - } - if(verboseconvolution)if(VerboseConvergence())_printf_("\n"); - - /*Free ressources:*/ - xDelete(Sg_old); -} -/*}}}*/ -void FemModel::SealevelriseRotationalFeedback(Vector* pSgo_rot, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/ - - /*serialized vectors:*/ - IssmDouble* Sg_old=NULL; - IssmDouble eartharea=0; - IssmDouble eartharea_cpu=0; - IssmDouble tide_love_h, tide_love_k, fluid_love, moi_e, moi_p, omega, g; - IssmDouble load_love_k2 = -0.30922675; //degree 2 load Love number - IssmDouble m1, m2, m3; - IssmDouble lati, longi, radi, value; - - /*Serialize vectors from previous iteration:*/ - Sg_old=pSg_old->ToMPISerial(); - - /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */ - for(int i=0;iSize();i++){ - 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()); - - IssmDouble moi_list[3]={0,0,0}; - IssmDouble moi_list_cpu[3]={0,0,0}; - for(int i=0;iSize();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->SealevelriseMomentOfInertia(&moi_list[0],Sg_old,eartharea); - moi_list_cpu[0] += moi_list[0]; - moi_list_cpu[1] += moi_list[1]; - moi_list_cpu[2] += moi_list[2]; - } - ISSM_MPI_Reduce (&moi_list_cpu[0],&moi_list[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&moi_list[0],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - // - ISSM_MPI_Reduce (&moi_list_cpu[1],&moi_list[1],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&moi_list[1],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - // - ISSM_MPI_Reduce (&moi_list_cpu[2],&moi_list[2],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&moi_list[2],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - - /*pull out some useful parameters: */ - parameters->FindParam(&tide_love_h,SealevelriseTidalLoveHEnum); - parameters->FindParam(&tide_love_k,SealevelriseTidalLoveKEnum); - parameters->FindParam(&fluid_love,SealevelriseFluidLoveEnum); - parameters->FindParam(&moi_e,SealevelriseEquatorialMoiEnum); - parameters->FindParam(&moi_p,SealevelrisePolarMoiEnum); - parameters->FindParam(&omega,SealevelriseAngularVelocityEnum); - - /*compute perturbation terms for angular velocity vector: */ - m1 = 1/(1-tide_love_k/fluid_love) * (1+load_love_k2)/(moi_p-moi_e) * moi_list[0]; - m2 = 1/(1-tide_love_k/fluid_love) * (1+load_love_k2)/(moi_p-moi_e) * moi_list[1]; - m3 = -(1+load_love_k2)/moi_p * moi_list[2]; // term associated with fluid number (3-order-of-magnitude smaller) is negelected - - /* Green's function (1+k_2-h_2/g): checked against Glenn Milne's thesis Chapter 3 (eqs: 3.3-4, 3.10-11) - * Perturbation terms for angular velocity vector (m1, m2, m3): checked against Mitrovica (2005 Appendix) & Jensen et al (2013 Appendix A3) - * Sea level rotational feedback: checked against GMD eqs 8-9 (only first order terms, i.e., degree 2 order 0 & 1 considered) - * all DONE in Geographic coordinates: theta \in [-90,90], lambda \in [-180 180] - */ - for(int i=0;iSize();i++){ - int sid; - //Vertex* vertex=(Vertex*)vertices->GetObjectByOffset(i); - Vertex* vertex=xDynamicCast(vertices->GetObjectByOffset(i)); - sid=vertex->Sid(); - - lati=latitude[sid]/180*PI; longi=longitude[sid]/180*PI; radi=radius[sid]; - - /*only first order terms are considered now: */ - value=((1.0+tide_love_k-tide_love_h)/9.81)*pow(omega*radi,2.0)* - (-m3/6.0 + 0.5*m3*cos(2.0*lati) - 0.5*sin(2.*lati)*(m1*cos(longi)+m2*sin(longi))); - - pSgo_rot->SetValue(sid,value,INS_VAL); //INS_VAL ensures that you don't add several times - } - - /*Assemble mesh velocity*/ - pSgo_rot->Assemble(); - - /*Free ressources:*/ - xDelete(Sg_old); - -} -/*}}}*/ -void FemModel::SealevelriseGeodetic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pSg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){/*{{{*/ - - /*serialized vectors:*/ - IssmDouble* Sg=NULL; - - IssmDouble eartharea=0; - IssmDouble eartharea_cpu=0; - - int ns,nsmax; - - /*Serialize vectors from previous iteration:*/ - Sg=pSg->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)); - 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()); - - /*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()); - - /*Call the sea level rise core: */ - for(int i=0;i(elements->GetObjectByOffset(i)); - element->SealevelriseGeodetic(pUp,pNorth,pEast,Sg,latitude,longitude,radius,xx,yy,zz,eartharea); - } - if(i%100==0){ - pUp->Assemble(); - pNorth->Assemble(); - pEast->Assemble(); - } - } - - /*One last time: */ - pUp->Assemble(); - pNorth->Assemble(); - pEast->Assemble(); - - /*Free ressources:*/ - xDelete(Sg); - xDelete(latitude); - xDelete(longitude); - xDelete(radius); - xDelete(xx); - xDelete(yy); - xDelete(zz); -} -/*}}}*/ -IssmDouble FemModel::SealevelriseOceanAverage(Vector* Sg) { /*{{{*/ - - IssmDouble* Sg_serial=NULL; - IssmDouble oceanvalue,oceanvalue_cpu; - IssmDouble oceanarea,oceanarea_cpu; - - /*Serialize vectors from previous iteration:*/ - Sg_serial=Sg->ToMPISerial(); - - /*Initialize:*/ - oceanvalue_cpu=0; - oceanarea_cpu=0; - - /*Go through elements, and add contribution from each element and divide by overall ocean area:*/ - for(int i=0;iSize();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - oceanarea_cpu += element->OceanArea(); - oceanvalue_cpu += element->OceanAverage(Sg_serial); - } - 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(Sg_serial); - - return oceanvalue/oceanarea; -} -/*}}}*/ -#endif -void FemModel::HydrologyEPLupdateDomainx(IssmDouble* pEplcount){ /*{{{*/ - - Vector* mask = NULL; - Vector* recurence = NULL; - Vector* active = NULL; - IssmDouble* serial_mask = NULL; - IssmDouble* serial_rec = NULL; - IssmDouble* serial_active = NULL; - IssmDouble* old_active = NULL; - int* eplzigzag_counter = NULL; - int eplflip_lock; - - HydrologyDCEfficientAnalysis* effanalysis = new HydrologyDCEfficientAnalysis(); - HydrologyDCInefficientAnalysis* inefanalysis = new HydrologyDCInefficientAnalysis(); - - /*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); - 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); - } - /*check for changes and increment zigzag counter, change the mask if necessary*/ - recurence->Assemble(); - serial_rec=recurence->ToMPISerial(); - for (int i=0;iSize();i++){ - Node* node=xDynamicCast(nodes->GetObjectByOffset(i)); - if(serial_rec[node->Sid()]==1.)eplzigzag_counter[node->Lid()] ++; - if(eplzigzag_counter[node->Lid()]>eplflip_lock & eplflip_lock!=0){ - mask->SetValue(node->Sid(),old_active[node->Sid()],INS_VAL); - } - } - this->parameters->SetParam(eplzigzag_counter,this->nodes->Size(),EplZigZagCounterEnum); - /*Assemble and serialize*/ - mask->Assemble(); - serial_mask=mask->ToMPISerial(); - - xDelete(eplzigzag_counter); - xDelete(serial_rec); - xDelete(old_active); - delete mask; - delete recurence; - - /*Update Mask*/ - InputUpdateFromVectorx(this,serial_mask,HydrologydcMaskEplactiveNodeEnum,NodeSIdEnum); - xDelete(serial_mask); - inefanalysis->ElementizeEplMask(this); - /*Step 2: update node activity. If one element is connected to mask=1, all nodes are active*/ - active=new Vector(nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum)); - for (int i=0;iSize();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - effanalysis->HydrologyEPLGetActive(active,element); - } - - /*Assemble and serialize*/ - active->Assemble(); - serial_active=active->ToMPISerial(); - delete active; - - /*Update node activation accordingly*/ - int counter =0; - for (int i=0;iSize();i++){ - Node* node=xDynamicCast(nodes->GetObjectByOffset(i)); - if(node->InAnalysis(HydrologyDCEfficientAnalysisEnum)){ - if(serial_active[node->Sid()]==1.){ - node->Activate(); - if(!node->IsClone()) counter++; - } - else{ - node->Deactivate(); - } - } - } - xDelete(serial_active); - delete effanalysis; - 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()); - counter=sum_counter; - *pEplcount = counter; - if(VerboseSolution()) _printf0_(" Number of active nodes in EPL layer: "<< counter <<"\n"); - - /*Update dof indexings*/ - this->UpdateConstraintsx(); - -} -/*}}}*/ -void FemModel::UpdateConstraintsL2ProjectionEPLx(IssmDouble* pL2count){ /*{{{*/ - - Vector* active = NULL; - IssmDouble* serial_active = NULL; - HydrologyDCEfficientAnalysis* effanalysis = new HydrologyDCEfficientAnalysis(); - - /*update node activity. If one element is connected to mask=1, all nodes are active*/ - active=new Vector(nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum)); - for (int i=0;iSize();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - effanalysis->HydrologyEPLGetActive(active,element); - } - - /*Assemble and serialize*/ - active->Assemble(); - serial_active=active->ToMPISerial(); - delete active; - delete effanalysis; - - /*Update node activation accordingly*/ - int counter =0; - for (int i=0;iSize();i++){ - Node* node=xDynamicCast(nodes->GetObjectByOffset(i)); - if(node->InAnalysis(L2ProjectionEPLAnalysisEnum)){ - if(serial_active[node->Sid()]==1.){ - node->Activate(); - if(!node->IsClone()) counter++; - } - else{ - node->Deactivate(); - } - } - } - 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()); - counter=sum_counter; - *pL2count = counter; - if(VerboseSolution()) _printf0_(" Number of active nodes L2 Projection: "<< counter <<"\n"); -} -/*}}}*/ -#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; - int ierr; - - /*First things first, store the communicator, and set it as a global variable: */ - IssmComm::SetComm(incomm); - - /*Start profiler: */ - this->profiler=new Profiler(); - profiler->Tag(START); - - /*From command line arguments, retrieve different filenames needed to create the FemModel: */ - solution_type=StringToEnumx(solution); - - /*Create femmodel from input files: */ - profiler->Tag(STARTINIT); - this->InitFromBuffers((char*)buffer,buffersize,toolkits, solution_type,trace,NULL); - profiler->Tag(FINISHINIT); - - /*Save communicator in the parameters dataset: */ - this->parameters->AddObject(new GenericParam(incomm,FemModelCommEnum)); - -} -/*}}}*/ -void FemModel::CleanUpJs(char** poutput, size_t* psize){/*{{{*/ - - /*Intermediary*/ - FILE *output_fid; - GenericParam* outputbufferparam=NULL; - GenericParam* outputbuffersizeparam=NULL; - char** poutputbuffer; - size_t* poutputbuffersize; - - - /*Before we delete the profiler, report statistics for this run: */ - profiler->Tag(FINISH); //final tagging - _printf0_("\n"); - _printf0_(" "<DeltaTime(STARTINIT,FINISHINIT) << "\n"); - _printf0_(" "<DeltaTime(STARTCORE,FINISHCORE) << "\n"); - _printf0_("\n"); - _printf0_(" Total elapsed time: " - <DeltaTimeModHour(START,FINISH)<<" hrs " - <DeltaTimeModMin(START,FINISH)<<" min " - <DeltaTimeModSec(START,FINISH)<<" sec" - ); - _printf0_("\n"); - - /*Before we close the output file, recover the buffer and size:*/ - outputbufferparam = xDynamicCast*>(this->parameters->FindParamObject(OutputBufferPointerEnum)); - poutputbuffer=outputbufferparam->GetParameterValue(); - outputbuffersizeparam = xDynamicCast*>(this->parameters->FindParamObject(OutputBufferSizePointerEnum)); - poutputbuffersize=outputbuffersizeparam->GetParameterValue(); - - /*Assign output values: */ - *poutput=*poutputbuffer; - *psize=*poutputbuffersize; -} -/*}}}*/ -void FemModel::InitFromBuffers(char* buffer, int buffersize, char* toolkits, int in_solution_type, bool trace, IssmPDouble* X){/*{{{*/ - - /*intermediary*/ - FILE *IOMODEL = NULL; - FILE *toolkitsoptionsfid = NULL; - FILE *output_fid = NULL; - int my_rank; - size_t outputsize; - char *outputbuffer; - const char *rootpath = ""; //needed for Dakota runs only, which we won't do here. - - /*recover my_rank:*/ - my_rank=IssmComm::GetRank(); - - /*Open input file descriptor on cpu 0: */ - if(my_rank==0) IOMODEL = fmemopen((void*)buffer, buffersize, "rb"); - - /*Open toolkits file descriptor: */ - toolkitsoptionsfid=fmemopen((void*)toolkits, strlen(toolkits)+1, "r"); - - /*Now, go create FemModel:*/ - this->InitFromFids((char*)rootpath,IOMODEL,toolkitsoptionsfid,in_solution_type,trace,X); - - /*Close input file and toolkits file descriptors: */ - if(my_rank==0) fclose(IOMODEL); - 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)_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)); - -}/*}}}*/ -#endif - -#ifdef _HAVE_NEOPZ_ -void FemModel::InitializeAdaptiveRefinement(void){/*{{{*/ - - /*Define variables*/ - int my_rank = IssmComm::GetRank(); - this->amr = NULL;//initialize amr as NULL - int numberofvertices = this->vertices->NumberOfVertices(); - int numberofelements = this->elements->NumberOfElements(); - int numberofsegments = 0; //used on matlab - IssmDouble* x = NULL; - IssmDouble* y = NULL; - IssmDouble* z = NULL; - int* elements = NULL; - int elementswidth = this->GetElementsWidth(); //just tria elements in this version. Itapopo: - int levelmax = 0; - IssmDouble regionlevel1 = 0.; - IssmDouble regionlevelmax = 0.; - - /*Get vertices coordinates of the coarse mesh (father mesh)*/ - /*elements comes in Matlab indexing*/ - this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elements); - - /*Get amr parameters*/ - this->parameters->FindParam(&levelmax,AmrLevelMaxEnum); - this->parameters->FindParam(®ionlevel1,AmrRegionLevel1Enum); - this->parameters->FindParam(®ionlevelmax,AmrRegionLevelMaxEnum); - - /*Create initial mesh (coarse mesh) in neopz data structure*/ - /*Just CPU #0 should keep AMR object*/ - this->SetRefPatterns(); - if(my_rank==0){ - bool ismisomip = false; - if(ismisomip){//itapopo - TPZFileStream fstr; - std::stringstream ss; - int frictionlaw; - - this->parameters->FindParam(&frictionlaw,FrictionLawEnum); - - ss << levelmax; - if(frictionlaw==1){ - ss << "_viscous/amr.txt"; - }else if(frictionlaw==7){ - ss << "_tsai/amr.txt"; - }else{ - _error_("friction law not supported here."); - } - - std::string AMRfile = "/home/santos/Misomip2/L" + ss.str(); - fstr.OpenRead(AMRfile.c_str()); - - TPZSaveable *sv = TPZSaveable::Restore(fstr,0); - this->amr = dynamic_cast(sv); - } - else{ - this->amr = new AdaptiveMeshRefinement(); - //this->amr->SetLevelMax(levelmax); //Set max level of refinement - //this->amr->SetRegions(regionlevel1,regionlevelmax); - this->amr->CreateInitialMesh(numberofvertices,numberofelements,numberofsegments,elementswidth,x,y,z,elements,NULL); - } - this->amr->SetLevelMax(levelmax); //Set max level of refinement - this->amr->SetRegions(regionlevel1,regionlevelmax); - } - - /*Free the vectors*/ - xDelete(x); - xDelete(y); - xDelete(z); - xDelete(elements); - -} -/*}}}*/ -void FemModel::SetRefPatterns(){/*{{{*/ - - /*Initialize the global variable of refinement patterns*/ - gRefDBase.InitializeUniformRefPattern(EOned); - gRefDBase.InitializeUniformRefPattern(ETriangle); - - //gRefDBase.InitializeRefPatterns(); - /*Insert specifics patterns to ISSM core*/ - std::string filepath = REFPATTERNDIR; - std::string filename1 = filepath + "/2D_Triang_Rib_3.rpt"; - std::string filename2 = filepath + "/2D_Triang_Rib_4.rpt"; - std::string filename3 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4.rpt"; - std::string filename4 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4_permuted.rpt"; - std::string filename5 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5.rpt"; - std::string filename6 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5_permuted.rpt"; - std::string filename7 = filepath + "/2D_Triang_Rib_5.rpt"; - - TPZAutoPointer refpat1 = new TPZRefPattern(filename1); - TPZAutoPointer refpat2 = new TPZRefPattern(filename2); - TPZAutoPointer refpat3 = new TPZRefPattern(filename3); - TPZAutoPointer refpat4 = new TPZRefPattern(filename4); - TPZAutoPointer refpat5 = new TPZRefPattern(filename5); - TPZAutoPointer refpat6 = new TPZRefPattern(filename6); - TPZAutoPointer refpat7 = new TPZRefPattern(filename7); - - if(!gRefDBase.FindRefPattern(refpat1)) gRefDBase.InsertRefPattern(refpat1); - if(!gRefDBase.FindRefPattern(refpat2)) gRefDBase.InsertRefPattern(refpat2); - if(!gRefDBase.FindRefPattern(refpat3)) gRefDBase.InsertRefPattern(refpat3); - if(!gRefDBase.FindRefPattern(refpat4)) gRefDBase.InsertRefPattern(refpat4); - if(!gRefDBase.FindRefPattern(refpat5)) gRefDBase.InsertRefPattern(refpat5); - if(!gRefDBase.FindRefPattern(refpat6)) gRefDBase.InsertRefPattern(refpat6); - if(!gRefDBase.FindRefPattern(refpat7)) gRefDBase.InsertRefPattern(refpat7); -} -/*}}}*/ -void FemModel::ReMesh(void){/*{{{*/ - - /*Variables*/ - IssmDouble *newx = NULL; - IssmDouble *newy = NULL; - IssmDouble *newz = NULL; - int *newelementslist = NULL; - int newnumberofvertices = -1; - int newnumberofelements = -1; - bool* my_elements = NULL; - int* my_vertices = NULL; - int elementswidth = this->GetElementsWidth();//just tria elements in this version - - /*Execute refinement and get the new mesh*/ - /*newelementslist come in Matlab indexing*/ - this->ExecuteRefinement(newnumberofvertices,newnumberofelements,&newx,&newy,&newz,&newelementslist); - - /*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); - - if(this->loads->Size()!=0) _error_("not supported yet"); - - /*Create vertices*/ - Vertices* new_vertices=new Vertices(); - this->CreateVertices(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,my_vertices,newx,newy,newz,new_vertices); - - /*Creating elements*/ - /*Just Tria in this version*/ - Elements* new_elements=new Elements(); - this->CreateElements(newnumberofelements,elementswidth,newelementslist,my_elements,new_elements); - - /*Creating materials*/ - Materials* new_materials=new Materials(); - this->CreateMaterials(newnumberofelements,my_elements,new_materials); - - /*Creating nodes and constraints*/ - /*Just SSA (2D) and P1 in this version*/ - Nodes* new_nodes=new Nodes(); - Constraints* new_constraints=new Constraints(); - - int nodecounter=0; - int constraintcounter=0; - for(int i=0;inummodels;i++){//create nodes for each analysis in analysis_type_list - - int analysis_enum = this->analysis_type_list[i]; - - /*As the domain is 2D, it is not necessary to create nodes for this analysis*/ - /*itapopo must verify if domain is not 3D. Only 2D in this version!*/ - if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue; - - this->CreateNodes(newnumberofvertices,my_vertices,nodecounter,analysis_enum,new_nodes); - if(analysis_enum==StressbalanceAnalysisEnum) this->CreateConstraints(newnumberofvertices,newnumberofelements,nodecounter,constraintcounter,newx,newy,my_vertices,new_constraints); - this->UpdateElements(newnumberofelements,newelementslist,my_elements,nodecounter,i,new_elements); - - if(new_nodes->Size()) nodecounter=new_nodes->MaximumId(); - constraintcounter = new_constraints->NumberOfConstraints(); - /*Make sure nodecounter is at least 0 (if no node exists, maxid will be -1*/ - _assert_(nodecounter>=0); - } - - new_elements->Presort(); - new_nodes->Presort(); - new_vertices->Presort(); - this->loads->Presort(); - new_materials->Presort(); - new_constraints->Presort(); - - /*reset hooks for elements, loads and nodes: */ - new_elements->ResetHooks(); - this->loads->ResetHooks(); - new_materials->ResetHooks(); - - /*do the post-processing of the datasets to get an FemModel that can actually run analyses: */ - int analysis_type; - for(int i=0;inummodels;i++){ - analysis_type=this->analysis_type_list[i]; - //SetCurrentConfiguration(analysis_type); - - 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); - this->parameters->SetParam(analysis_type,ConfigurationTypeEnum); - - /*configure elements, loads and nodes, for this new analysis: */ - new_elements->SetCurrentConfiguration(new_elements,this->loads,new_nodes,new_vertices,new_materials,this->parameters); - this->loads->SetCurrentConfiguration(new_elements,this->loads,new_nodes,new_vertices,new_materials,this->parameters); - - /*take care of toolkits options, that depend on this analysis type (present only after model processor)*/ - if(this->parameters->Exist(ToolkitsOptionsStringsEnum)){ - ToolkitsOptionsFromAnalysis(this->parameters,analysis_type); - if(VerboseSolver()) _printf0_(" toolkits Options set for analysis type: " << EnumToStringx(analysis_type) << "\n"); - } - - ConfigureObjectsx(new_elements,this->loads,new_nodes,new_vertices,new_materials,this->parameters); - 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); - NodesDofx(new_nodes,this->parameters,analysis_type); - } - - /*Finally: interpolate all inputs and insert them into the new elements.*/ - this->InterpolateInputs(new_vertices,new_elements); - - /*Delete old structure and set new pointers*/ - delete this->vertices; this->vertices = new_vertices; - delete this->elements; this->elements = new_elements; - delete this->nodes; this->nodes = new_nodes; - delete this->constraints; this->constraints = new_constraints; - delete this->materials; this->materials = new_materials; - - GetMaskOfIceVerticesLSMx(this); - - /*Insert MISMIP+ bed topography*/ - if(false) this->BedrockFromMismipPlus(); - - /*Adjust base, thickness and mask grounded ice leve set*/ - if(true) this->AdjustBaseThicknessAndMask(); - - /*Reset current configuration: */ - analysis_type=this->analysis_type_list[this->analysis_counter]; - SetCurrentConfiguration(analysis_type); - - /*Cleanup*/ - xDelete(newx); - xDelete(newy); - xDelete(newz); - xDelete(newelementslist); - xDelete(my_vertices); - xDelete(my_elements); - - return; - -} -/*}}}*/ + +/*AMR*/ void FemModel::BedrockFromMismipPlus(void){/*{{{*/ @@ -3245,5 +2466,5 @@ } - if(abs(sl[i])>0) _error_("Sea level value not supported!"); + if(abs(sl[i])>0) _error_("Sea level value "<AddInput(ThicknessEnum,&h[0],P1Enum); element->AddInput(BaseEnum,&b[0],P1Enum); - } @@ -3265,68 +2485,4 @@ xDelete(r); xDelete(sl); - - return; -} -/*}}}*/ -void FemModel::WriteMeshInResults(void){/*{{{*/ - - int step = -1; - int numberofelements = -1; - int numberofvertices = -1; - IssmDouble time = -1; - IssmDouble* x = NULL; - IssmDouble* y = NULL; - IssmDouble* z = NULL; - 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(); - numberofvertices=this->vertices->NumberOfVertices(); - - /*Get mesh. Elementslist comes in Matlab indexing*/ - this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elementslist); - - /*Write mesh in Results*/ - this->results->AddResult(new GenericExternalResult(this->results->Size()+1,MeshElementsEnum, - elementslist,numberofelements,this->GetElementsWidth(),step,time)); - - this->results->AddResult(new GenericExternalResult(this->results->Size()+1,MeshXEnum, - x,numberofvertices,1,step,time)); - - this->results->AddResult(new GenericExternalResult(this->results->Size()+1,MeshYEnum, - y,numberofvertices,1,step,time)); - - //itapopo - if(IssmComm::GetRank()==0){ - TPZFileStream fstr; - std::stringstream ss; - int frictionlaw,levelmax; - this->parameters->FindParam(&frictionlaw,FrictionLawEnum); - this->parameters->FindParam(&levelmax,AmrLevelMaxEnum); - ss << levelmax; - if(frictionlaw==1){ - ss << "_viscous/amr.txt"; - }else if(frictionlaw==7){ - ss << "_tsai/amr.txt"; - }else{ - _error_("friction law not supported here."); - } - std::string AMRfile = "/home/santos/Misomip2/L" + ss.str(); - fstr.OpenWrite(AMRfile.c_str()); - int withclassid = 1; - this->amr->Write(fstr,withclassid); - } - //itapopo - - /*Cleanup*/ - xDelete(x); - xDelete(y); - xDelete(z); - xDelete(elementslist); - - return; } /*}}}*/ @@ -3342,5 +2498,5 @@ } input_interpolations->Assemble(); - + /*Serialize and set output*/ IssmDouble* input_interpolations_serial = input_interpolations->ToMPISerial(); @@ -3374,5 +2530,5 @@ P1input_interp[numP1inputs] = inputinterp; } - numP1inputs++; + numP1inputs++; break; case P0Enum: @@ -3421,5 +2577,5 @@ xDelete(vector); } - + /*Old mesh coordinates*/ IssmDouble *Xold = NULL; @@ -3433,5 +2589,5 @@ IssmDouble* YC_new = NULL; int *Indexnew = NULL; - + /*Get the old mesh*/ this->GetMesh(this->vertices,this->elements,&Xold,&Yold,&Zold,&Indexold); @@ -3459,14 +2615,14 @@ P0inputsold,numelementsold,numP0inputs, XC_new,YC_new,numelementsnew,NULL); - + /*Interpolate P1 inputs in the new mesh*/ InterpFromMeshToMesh2dx(&P1inputsnew,Indexold,Xold,Yold,numverticesold,numelementsold, P1inputsold,numverticesold,numP1inputs, Xnew,Ynew,numverticesnew,NULL); - + /*Insert P0 inputs into the new elements.*/ vector=NULL; for(int i=0;i(numelementsnew); @@ -3533,65 +2689,67 @@ xDelete(YC_new); xDelete(Indexnew); - -} -/*}}}*/ -void FemModel::ExecuteRefinement(int &numberofvertices,int &numberofelements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist){/*{{{*/ - - /*elements is in Matlab indexing*/ - int my_rank = IssmComm::GetRank(); - int numberofsegments = -1; - IssmDouble* vx = NULL; //itapopo this is not being used - IssmDouble* vy = NULL; //itapopo this is not being used - IssmDouble* x = NULL; - IssmDouble* y = NULL; - IssmDouble* z = NULL; - int* elementslist = NULL; - int* segments = NULL; - IssmDouble* masklevelset = NULL; - IssmDouble* pelementerror = NULL; - const int elementswidth = this->GetElementsWidth();//just 2D mesh, tria elements - - /*Solutions which will be used to refine the elements*/ - this->GetGroundediceLevelSet(&masklevelset);//itapopo verificar se já existe um método igual a esse - - //Compute the ZZ error estimator per element - this->ZZErrorEstimator(&pelementerror); - - _printf0_("P Element error\n"); - for(int i=0;ielements->NumberOfElements();i++) _printf0_(""<amr->ExecuteRefinement(type_process,vx,vy,masklevelset, - numberofvertices,numberofelements,numberofsegments,&x,&y,&z,&elementslist,&segments); - if(numberofvertices<=0 || numberofelements<=0 /*|| newnumberofsegments<=0*/) _error_("Error in the refinement process."); - } - else{ - x=xNew(numberofvertices); - y=xNew(numberofvertices); - z=xNew(numberofvertices); - elementslist=xNew(numberofelements*this->GetElementsWidth()); - } - - /*Send new mesh to others CPU*/ - ISSM_MPI_Bcast(&numberofvertices,1,ISSM_MPI_INT,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(&numberofelements,1,ISSM_MPI_INT,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(x,numberofvertices,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(y,numberofvertices,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(z,numberofvertices,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(elementslist,numberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm()); - - /*Assign the pointers*/ - (*pelementslist) = elementslist; //Matlab indexing - (*px) = x; - (*py) = y; - (*pz) = z; - +} +/*}}}*/ +void FemModel::WriteMeshInResults(void){/*{{{*/ + + int step = -1; + int numberofelements = -1; + int numberofvertices = -1; + IssmDouble time = -1; + IssmDouble* x = NULL; + IssmDouble* y = NULL; + IssmDouble* z = NULL; + 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(); + numberofvertices=this->vertices->NumberOfVertices(); + + /*Get mesh. Elementslist comes in Matlab indexing*/ + this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elementslist); + + /*Write mesh in Results*/ + printf("-------------- file: FemModel.cpp line: %i\n",__LINE__); + this->results->AddResult(new GenericExternalResult(this->results->Size()+1,MeshElementsEnum, + elementslist,numberofelements,this->GetElementsWidth(),step,time)); + + this->results->AddResult(new GenericExternalResult(this->results->Size()+1,MeshXEnum, + x,numberofvertices,1,step,time)); + + this->results->AddResult(new GenericExternalResult(this->results->Size()+1,MeshYEnum, + y,numberofvertices,1,step,time)); + + //itapopo + #ifdef _HAVE_NEOPZ_ + if(IssmComm::GetRank()==0){ + TPZFileStream fstr; + std::stringstream ss; + int frictionlaw,levelmax; + this->parameters->FindParam(&frictionlaw,FrictionLawEnum); + this->parameters->FindParam(&levelmax,AmrLevelMaxEnum); + ss << levelmax; + if(frictionlaw==1){ + ss << "_viscous/amr.txt"; + }else if(frictionlaw==7){ + ss << "_tsai/amr.txt"; + }else{ + _error_("friction law not supported here."); + } + std::string AMRfile = "/home/santos/Misomip2/L" + ss.str(); + fstr.OpenWrite(AMRfile.c_str()); + int withclassid = 1; + this->amr->Write(fstr,withclassid); + } + #endif + //itapopo + /*Cleanup*/ - if(segments) xDelete(segments); - xDelete(masklevelset); - xDelete(pelementerror); - + xDelete(x); + xDelete(y); + xDelete(z); + xDelete(elementslist); } /*}}}*/ @@ -3635,6 +2793,6 @@ connectivity=xNewZeroInit(newnumberofvertices); - for (int i=0;i0 && vertexid-1SetMid(newnumberofelements+1); materials->AddObject(newmatpar);//put it at the end of the materials - } /*}}}*/ @@ -3954,5 +3111,4 @@ delete vspcvxflag; delete vspcvyflag; - } /*}}}*/ @@ -4053,5 +3209,4 @@ xDelete(npart); xDelete(index); - } /*}}}*/ @@ -4141,5 +3296,4 @@ xDelete(totalweight); xDelete(elem_vertices); - } /*}}}*/ @@ -4215,6 +3369,964 @@ xDelete(elem_vertices); delete velementerror; +} +/*}}}*/ + +#ifdef _HAVE_DAKOTA_ +void FemModel::DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses){/*{{{*/ + + int i,j; + int my_rank; + + /*intermediary: */ + char root[50]; + int index; + int npart; + double femmodel_response; + int flag; + double *vertex_response = NULL; + double *qmu_response = NULL; + double *responses_pointer = NULL; + + /*retrieve npart: */ + parameters->FindParam(&npart,QmuNumberofpartitionsEnum); + my_rank=IssmComm::GetRank(); + + /*save the d_responses pointer: */ + responses_pointer=d_responses; + + //watch out, we have more d_numresponses than numresponsedescriptors, because the responses have been expanded if they were scaled. + //because we don't know the d_responses descriptors (the scaled ones) we can't key off them, so we will key off the responses_descriptors: */ + + for(i=0;i(vertex_response); + xDelete(qmu_response); + + } + else if (flag==IndexedEnum){ + + /*indexed response: plug index into parameters and call response module: */ + parameters->SetParam(index,IndexEnum); + + this->Responsex(&femmodel_response,root); + + if(my_rank==0){ + /*plug response: */ + responses_pointer[0]=femmodel_response; + + /*increment response_pointer :*/ + responses_pointer++; + } + } + else if (flag==NodalEnum){ + _error_("nodal response functions not supported yet!"); + + /*increment response_pointer :*/ + responses_pointer++; + } + else if (flag==RegularEnum){ + + /*perfectly normal response function: */ + this->Responsex(&femmodel_response,root); + + if(my_rank==0){ + /*plug response: */ + responses_pointer[0]=femmodel_response; + + /*increment response_pointer :*/ + responses_pointer++; + } + } + else _error_("flag type " << flag << " not supported yet for response analysis"); + } + + /*Synthesize echo: {{{*/ + if(my_rank==0){ + _printf_(" responses: " << d_numresponses << ": "); + for(i=0;i* wg,Vector* dwgdt, IssmDouble* x, IssmDouble* y){ /*{{{*/ + + /*Go through elements, and add contribution from each element to the deflection vector wg:*/ + for(int i=0;iSize();i++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + element->GiaDeflection(wg,dwgdt, x,y); + } +} +/*}}}*/ +#endif +#ifdef _HAVE_ESA_ +void FemModel::EsaGeodetic2D(Vector* pUp, Vector* pNorth, Vector* pEast, 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()); + + /*Call the esa geodetic core: */ + for(int i=0;i(elements->GetObjectByOffset(i)); + element->EsaGeodetic2D(pUp,pNorth,pEast,xx,yy); + } + if(i%100==0){ + pUp->Assemble(); + pNorth->Assemble(); + pEast->Assemble(); + } + } + + /*One last time: */ + pUp->Assemble(); + pNorth->Assemble(); + pEast->Assemble(); + + /*Free ressources:*/ + xDelete(xx); + xDelete(yy); +} +/*}}}*/ +void FemModel::EsaGeodetic3D(Vector* pUp, Vector* pNorth, Vector* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){/*{{{*/ + + IssmDouble eartharea=0; + 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: */ + for(int i=0;i(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()); + + /*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()); + + /*Call the esa geodetic core: */ + for(int i=0;i(elements->GetObjectByOffset(i)); + element->EsaGeodetic3D(pUp,pNorth,pEast,latitude,longitude,radius,xx,yy,zz,eartharea); + } + if(i%100==0){ + pUp->Assemble(); + pNorth->Assemble(); + pEast->Assemble(); + } + } + + /*One last time: */ + pUp->Assemble(); + pNorth->Assemble(); + pEast->Assemble(); + + /*Free ressources:*/ + xDelete(latitude); + xDelete(longitude); + xDelete(radius); + xDelete(xx); + xDelete(yy); + xDelete(zz); +} +/*}}}*/ +#endif +#ifdef _HAVE_SEALEVELRISE_ +void FemModel::SealevelriseEustatic(Vector* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius) { /*{{{*/ + + /*serialized vectors:*/ + IssmDouble eustatic = 0.; + IssmDouble eustatic_cpu = 0.; + IssmDouble eustatic_cpu_e = 0.; + IssmDouble oceanarea = 0.; + IssmDouble oceanarea_cpu = 0.; + IssmDouble eartharea = 0.; + 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 area of the ocean, which is needed to compute the eustatic component: */ + for(int i=0;i(elements->GetObjectByOffset(i)); + oceanarea_cpu += element->OceanArea(); + eartharea_cpu += element->GetAreaSpherical(); + } + 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()); + _assert_(oceanarea>0.); + + 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()); + + /*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()); + + /*Call the sea level rise core: */ + for(int i=0;i(elements->GetObjectByOffset(i)); + element->SealevelriseEustatic(pSgi,&eustatic_cpu_e,latitude,longitude,radius,oceanarea,eartharea); + eustatic_cpu+=eustatic_cpu_e; + } + if(i%100==0)pSgi->Assemble(); + } + if(VerboseConvergence())_printf0_("\n"); + + /*One last time: */ + pSgi->Assemble(); + + /*Sum all eustatic components from all cpus:*/ + ISSM_MPI_Reduce (&eustatic_cpu,&eustatic,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&eustatic,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + _assert_(!xIsNan(eustatic)); + + /*Assign output pointers:*/ + *peustatic=eustatic; + +} +/*}}}*/ +void FemModel::SealevelriseNonEustatic(Vector* pSgo, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution){/*{{{*/ + + /*serialized vectors:*/ + IssmDouble* Sg_old=NULL; + + IssmDouble eartharea=0; + IssmDouble eartharea_cpu=0; + + int ns,nsmax; + + /*Serialize vectors from previous iteration:*/ + Sg_old=pSg_old->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)); + 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()); + + /*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()); + + /*Call the sea level rise core: */ + for(int i=0;i(elements->GetObjectByOffset(i)); + element->SealevelriseNonEustatic(pSgo,Sg_old,latitude,longitude,radius,eartharea); + } + if(i%100==0)pSgo->Assemble(); + } + if(verboseconvolution)if(VerboseConvergence())_printf_("\n"); + + /*Free ressources:*/ + xDelete(Sg_old); +} +/*}}}*/ +void FemModel::SealevelriseRotationalFeedback(Vector* pSgo_rot, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/ + + /*serialized vectors:*/ + IssmDouble* Sg_old=NULL; + IssmDouble eartharea=0; + IssmDouble eartharea_cpu=0; + IssmDouble tide_love_h, tide_love_k, fluid_love, moi_e, moi_p, omega, g; + IssmDouble load_love_k2 = -0.30922675; //degree 2 load Love number + IssmDouble m1, m2, m3; + IssmDouble lati, longi, radi, value; + + /*Serialize vectors from previous iteration:*/ + Sg_old=pSg_old->ToMPISerial(); + + /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */ + for(int i=0;iSize();i++){ + 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()); + + IssmDouble moi_list[3]={0,0,0}; + IssmDouble moi_list_cpu[3]={0,0,0}; + for(int i=0;iSize();i++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + element->SealevelriseMomentOfInertia(&moi_list[0],Sg_old,eartharea); + moi_list_cpu[0] += moi_list[0]; + moi_list_cpu[1] += moi_list[1]; + moi_list_cpu[2] += moi_list[2]; + } + ISSM_MPI_Reduce (&moi_list_cpu[0],&moi_list[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&moi_list[0],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + // + ISSM_MPI_Reduce (&moi_list_cpu[1],&moi_list[1],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&moi_list[1],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + // + ISSM_MPI_Reduce (&moi_list_cpu[2],&moi_list[2],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&moi_list[2],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + + /*pull out some useful parameters: */ + parameters->FindParam(&tide_love_h,SealevelriseTidalLoveHEnum); + parameters->FindParam(&tide_love_k,SealevelriseTidalLoveKEnum); + parameters->FindParam(&fluid_love,SealevelriseFluidLoveEnum); + parameters->FindParam(&moi_e,SealevelriseEquatorialMoiEnum); + parameters->FindParam(&moi_p,SealevelrisePolarMoiEnum); + parameters->FindParam(&omega,SealevelriseAngularVelocityEnum); + + /*compute perturbation terms for angular velocity vector: */ + m1 = 1/(1-tide_love_k/fluid_love) * (1+load_love_k2)/(moi_p-moi_e) * moi_list[0]; + m2 = 1/(1-tide_love_k/fluid_love) * (1+load_love_k2)/(moi_p-moi_e) * moi_list[1]; + m3 = -(1+load_love_k2)/moi_p * moi_list[2]; // term associated with fluid number (3-order-of-magnitude smaller) is negelected + + /* Green's function (1+k_2-h_2/g): checked against Glenn Milne's thesis Chapter 3 (eqs: 3.3-4, 3.10-11) + * Perturbation terms for angular velocity vector (m1, m2, m3): checked against Mitrovica (2005 Appendix) & Jensen et al (2013 Appendix A3) + * Sea level rotational feedback: checked against GMD eqs 8-9 (only first order terms, i.e., degree 2 order 0 & 1 considered) + * all DONE in Geographic coordinates: theta \in [-90,90], lambda \in [-180 180] + */ + for(int i=0;iSize();i++){ + int sid; + //Vertex* vertex=(Vertex*)vertices->GetObjectByOffset(i); + Vertex* vertex=xDynamicCast(vertices->GetObjectByOffset(i)); + sid=vertex->Sid(); + + lati=latitude[sid]/180*PI; longi=longitude[sid]/180*PI; radi=radius[sid]; + + /*only first order terms are considered now: */ + value=((1.0+tide_love_k-tide_love_h)/9.81)*pow(omega*radi,2.0)* + (-m3/6.0 + 0.5*m3*cos(2.0*lati) - 0.5*sin(2.*lati)*(m1*cos(longi)+m2*sin(longi))); + + pSgo_rot->SetValue(sid,value,INS_VAL); //INS_VAL ensures that you don't add several times + } + + /*Assemble mesh velocity*/ + pSgo_rot->Assemble(); + + /*Free ressources:*/ + xDelete(Sg_old); + +} +/*}}}*/ +void FemModel::SealevelriseGeodetic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pSg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){/*{{{*/ + + /*serialized vectors:*/ + IssmDouble* Sg=NULL; + + IssmDouble eartharea=0; + IssmDouble eartharea_cpu=0; + + int ns,nsmax; + + /*Serialize vectors from previous iteration:*/ + Sg=pSg->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)); + 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()); + + /*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()); + + /*Call the sea level rise core: */ + for(int i=0;i(elements->GetObjectByOffset(i)); + element->SealevelriseGeodetic(pUp,pNorth,pEast,Sg,latitude,longitude,radius,xx,yy,zz,eartharea); + } + if(i%100==0){ + pUp->Assemble(); + pNorth->Assemble(); + pEast->Assemble(); + } + } + + /*One last time: */ + pUp->Assemble(); + pNorth->Assemble(); + pEast->Assemble(); + + /*Free ressources:*/ + xDelete(Sg); + xDelete(latitude); + xDelete(longitude); + xDelete(radius); + xDelete(xx); + xDelete(yy); + xDelete(zz); +} +/*}}}*/ +IssmDouble FemModel::SealevelriseOceanAverage(Vector* Sg) { /*{{{*/ + + IssmDouble* Sg_serial=NULL; + IssmDouble oceanvalue,oceanvalue_cpu; + IssmDouble oceanarea,oceanarea_cpu; + + /*Serialize vectors from previous iteration:*/ + Sg_serial=Sg->ToMPISerial(); + + /*Initialize:*/ + oceanvalue_cpu=0; + oceanarea_cpu=0; + + /*Go through elements, and add contribution from each element and divide by overall ocean area:*/ + for(int i=0;iSize();i++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + oceanarea_cpu += element->OceanArea(); + oceanvalue_cpu += element->OceanAverage(Sg_serial); + } + 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(Sg_serial); + + return oceanvalue/oceanarea; +} +/*}}}*/ +#endif +void FemModel::HydrologyEPLupdateDomainx(IssmDouble* pEplcount){ /*{{{*/ + + Vector* mask = NULL; + Vector* recurence = NULL; + Vector* active = NULL; + IssmDouble* serial_mask = NULL; + IssmDouble* serial_rec = NULL; + IssmDouble* serial_active = NULL; + IssmDouble* old_active = NULL; + int* eplzigzag_counter = NULL; + int eplflip_lock; + + HydrologyDCEfficientAnalysis* effanalysis = new HydrologyDCEfficientAnalysis(); + HydrologyDCInefficientAnalysis* inefanalysis = new HydrologyDCInefficientAnalysis(); + + /*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); + 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); + } + /*check for changes and increment zigzag counter, change the mask if necessary*/ + recurence->Assemble(); + serial_rec=recurence->ToMPISerial(); + for (int i=0;iSize();i++){ + Node* node=xDynamicCast(nodes->GetObjectByOffset(i)); + if(serial_rec[node->Sid()]==1.)eplzigzag_counter[node->Lid()] ++; + if(eplzigzag_counter[node->Lid()]>eplflip_lock & eplflip_lock!=0){ + mask->SetValue(node->Sid(),old_active[node->Sid()],INS_VAL); + } + } + this->parameters->SetParam(eplzigzag_counter,this->nodes->Size(),EplZigZagCounterEnum); + /*Assemble and serialize*/ + mask->Assemble(); + serial_mask=mask->ToMPISerial(); + + xDelete(eplzigzag_counter); + xDelete(serial_rec); + xDelete(old_active); + delete mask; + delete recurence; + + /*Update Mask*/ + InputUpdateFromVectorx(this,serial_mask,HydrologydcMaskEplactiveNodeEnum,NodeSIdEnum); + xDelete(serial_mask); + inefanalysis->ElementizeEplMask(this); + /*Step 2: update node activity. If one element is connected to mask=1, all nodes are active*/ + active=new Vector(nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum)); + for (int i=0;iSize();i++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + effanalysis->HydrologyEPLGetActive(active,element); + } + + /*Assemble and serialize*/ + active->Assemble(); + serial_active=active->ToMPISerial(); + delete active; + + /*Update node activation accordingly*/ + int counter =0; + for (int i=0;iSize();i++){ + Node* node=xDynamicCast(nodes->GetObjectByOffset(i)); + if(node->InAnalysis(HydrologyDCEfficientAnalysisEnum)){ + if(serial_active[node->Sid()]==1.){ + node->Activate(); + if(!node->IsClone()) counter++; + } + else{ + node->Deactivate(); + } + } + } + xDelete(serial_active); + delete effanalysis; + 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()); + counter=sum_counter; + *pEplcount = counter; + if(VerboseSolution()) _printf0_(" Number of active nodes in EPL layer: "<< counter <<"\n"); + + /*Update dof indexings*/ + this->UpdateConstraintsx(); + +} +/*}}}*/ +void FemModel::UpdateConstraintsL2ProjectionEPLx(IssmDouble* pL2count){ /*{{{*/ + + Vector* active = NULL; + IssmDouble* serial_active = NULL; + HydrologyDCEfficientAnalysis* effanalysis = new HydrologyDCEfficientAnalysis(); + + /*update node activity. If one element is connected to mask=1, all nodes are active*/ + active=new Vector(nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum)); + for (int i=0;iSize();i++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + effanalysis->HydrologyEPLGetActive(active,element); + } + + /*Assemble and serialize*/ + active->Assemble(); + serial_active=active->ToMPISerial(); + delete active; + delete effanalysis; + + /*Update node activation accordingly*/ + int counter =0; + for (int i=0;iSize();i++){ + Node* node=xDynamicCast(nodes->GetObjectByOffset(i)); + if(node->InAnalysis(L2ProjectionEPLAnalysisEnum)){ + if(serial_active[node->Sid()]==1.){ + node->Activate(); + if(!node->IsClone()) counter++; + } + else{ + node->Deactivate(); + } + } + } + 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()); + counter=sum_counter; + *pL2count = counter; + if(VerboseSolution()) _printf0_(" Number of active nodes L2 Projection: "<< counter <<"\n"); +} +/*}}}*/ +#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; + int ierr; + + /*First things first, store the communicator, and set it as a global variable: */ + IssmComm::SetComm(incomm); + + /*Start profiler: */ + this->profiler=new Profiler(); + profiler->Tag(START); + + /*From command line arguments, retrieve different filenames needed to create the FemModel: */ + solution_type=StringToEnumx(solution); + + /*Create femmodel from input files: */ + profiler->Tag(STARTINIT); + this->InitFromBuffers((char*)buffer,buffersize,toolkits, solution_type,trace,NULL); + profiler->Tag(FINISHINIT); + + /*Save communicator in the parameters dataset: */ + this->parameters->AddObject(new GenericParam(incomm,FemModelCommEnum)); + +} +/*}}}*/ +void FemModel::CleanUpJs(char** poutput, size_t* psize){/*{{{*/ + + /*Intermediary*/ + FILE *output_fid; + GenericParam* outputbufferparam=NULL; + GenericParam* outputbuffersizeparam=NULL; + char** poutputbuffer; + size_t* poutputbuffersize; + + + /*Before we delete the profiler, report statistics for this run: */ + profiler->Tag(FINISH); //final tagging + _printf0_("\n"); + _printf0_(" "<DeltaTime(STARTINIT,FINISHINIT) << "\n"); + _printf0_(" "<DeltaTime(STARTCORE,FINISHCORE) << "\n"); + _printf0_("\n"); + _printf0_(" Total elapsed time: " + <DeltaTimeModHour(START,FINISH)<<" hrs " + <DeltaTimeModMin(START,FINISH)<<" min " + <DeltaTimeModSec(START,FINISH)<<" sec" + ); + _printf0_("\n"); + + /*Before we close the output file, recover the buffer and size:*/ + outputbufferparam = xDynamicCast*>(this->parameters->FindParamObject(OutputBufferPointerEnum)); + poutputbuffer=outputbufferparam->GetParameterValue(); + outputbuffersizeparam = xDynamicCast*>(this->parameters->FindParamObject(OutputBufferSizePointerEnum)); + poutputbuffersize=outputbuffersizeparam->GetParameterValue(); + + /*Assign output values: */ + *poutput=*poutputbuffer; + *psize=*poutputbuffersize; +} +/*}}}*/ +void FemModel::InitFromBuffers(char* buffer, int buffersize, char* toolkits, int in_solution_type, bool trace, IssmPDouble* X){/*{{{*/ + + /*intermediary*/ + FILE *IOMODEL = NULL; + FILE *toolkitsoptionsfid = NULL; + FILE *output_fid = NULL; + int my_rank; + size_t outputsize; + char *outputbuffer; + const char *rootpath = ""; //needed for Dakota runs only, which we won't do here. + + /*recover my_rank:*/ + my_rank=IssmComm::GetRank(); + + /*Open input file descriptor on cpu 0: */ + if(my_rank==0) IOMODEL = fmemopen((void*)buffer, buffersize, "rb"); + + /*Open toolkits file descriptor: */ + toolkitsoptionsfid=fmemopen((void*)toolkits, strlen(toolkits)+1, "r"); + + /*Now, go create FemModel:*/ + this->InitFromFids((char*)rootpath,IOMODEL,toolkitsoptionsfid,in_solution_type,trace,X); + + /*Close input file and toolkits file descriptors: */ + if(my_rank==0) fclose(IOMODEL); + 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)_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)); + +}/*}}}*/ +#endif + +#ifdef _HAVE_BAMG_ +void FemModel::ReMeshBamg(int* pnewnumberofvertices,int* pnewnumberofelements,IssmDouble** pnewx,IssmDouble** pnewy,IssmDouble** pnewz,int** pnewelementslist){/*{{{*/ + + /*Output*/ + IssmDouble *newx = NULL; + IssmDouble *newy = NULL; + IssmDouble *newz = NULL; + int *newelementslist = NULL; + int newnumberofvertices = -1; + int newnumberofelements = -1; + + /*Get Rank*/ + int my_rank = IssmComm::GetRank(); + + if(my_rank==0){ + this->amrbamg->ExecuteRefinementBamg(&newnumberofvertices,&newnumberofelements,&newx,&newy,&newz,&newelementslist); + if(newnumberofvertices<=0 || newnumberofelements<=0) _error_("Error in the refinement process."); + } + else{ + newx=xNew(newnumberofvertices); + newy=xNew(newnumberofvertices); + newz=xNew(newnumberofvertices); + newelementslist=xNew(newnumberofelements*this->GetElementsWidth()); + } + + /*Send new mesh to others CPU*/ + ISSM_MPI_Bcast(&newnumberofvertices,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&newnumberofelements,1,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; + *pnewnumberofelements = newnumberofelements; + *pnewx = newx; + *pnewy = newy; + *pnewz = newz; + *pnewelementslist = newelementslist; +} +/*}}}*/ +void FemModel::InitializeAdaptiveRefinementBamg(void){/*{{{*/ + + /*Define variables*/ + int numberofvertices = this->vertices->NumberOfVertices(); + int numberofelements = this->elements->NumberOfElements(); + int numberofsegments = 0; //used on matlab + IssmDouble* x = NULL; + IssmDouble* y = NULL; + IssmDouble* z = NULL; + int* elements = NULL; + int elementswidth = this->GetElementsWidth(); //just tria elements in this version. Itapopo: + IssmDouble hmin,hmax,err; + int fieldenum; + + /*Get rank*/ + int my_rank = IssmComm::GetRank(); + + /*Initialize field as NULL for now*/ + this->amrbamg = NULL; + + /*Get vertices coordinates of the coarse mesh (father mesh)*/ + this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elements); + + /*Get amr parameters*/ + this->parameters->FindParam(&hmin,AmrHminEnum); + this->parameters->FindParam(&hmax,AmrHmaxEnum); + this->parameters->FindParam(&fieldenum,AmrFieldEnum); + this->parameters->FindParam(&err,AmrErrEnum); + + /*Create bamg data structures for bamg (only cpu 0)*/ + if(my_rank==0){ + /*Re-create original mesh and put it in bamg structure*/ + this->amrbamg = new AmrBamg(hmin,hmax); + this->amrbamg->Initialize(elements,x,y,numberofvertices,numberofelements); + } + + /*Free the vectors*/ + xDelete(x); + xDelete(y); + xDelete(z); + xDelete(elements); +} +/*}}}*/ +#endif + +#ifdef _HAVE_NEOPZ_ +void FemModel::ReMeshNeopz(int* pnumberofvertices,int* pnumberofelements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist){/*{{{*/ + + /*elements is in Matlab indexing*/ + int my_rank = IssmComm::GetRank(); + int numberofsegments = -1; + int numberofvertices,numberofelements; + IssmDouble* vx = NULL; //itapopo this is not being used + IssmDouble* vy = NULL; //itapopo this is not being used + IssmDouble* x = NULL; + IssmDouble* y = NULL; + IssmDouble* z = NULL; + int* elementslist = NULL; + int* segments = NULL; + IssmDouble* masklevelset = NULL; + IssmDouble* pelementerror = NULL; + const int elementswidth = this->GetElementsWidth();//just 2D mesh, tria elements + + /*Solutions which will be used to refine the elements*/ + this->GetGroundediceLevelSet(&masklevelset);//itapopo verificar se já existe um método igual a esse + + //Compute the ZZ error estimator per element + this->ZZErrorEstimator(&pelementerror); + + _printf0_("P Element error\n"); + for(int i=0;ielements->NumberOfElements();i++) _printf0_(""<amr->ExecuteRefinement(type_process,vx,vy,masklevelset, + numberofvertices,numberofelements,numberofsegments,&x,&y,&z,&elementslist,&segments); + if(numberofvertices<=0 || numberofelements<=0 /*|| newnumberofsegments<=0*/) _error_("Error in the refinement process."); + } + else{ + x=xNew(numberofvertices); + y=xNew(numberofvertices); + z=xNew(numberofvertices); + elementslist=xNew(numberofelements*this->GetElementsWidth()); + } + + /*Send new mesh to others CPU*/ + ISSM_MPI_Bcast(&numberofvertices,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&numberofelements,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(x,numberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(y,numberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(z,numberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(elementslist,numberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm()); + + /*Assign the pointers*/ + (*pelementslist) = elementslist; //Matlab indexing + (*px) = x; + (*py) = y; + (*pz) = z; + *pnumberofelements = numberofelements; + *pnumberofvertices = numberofvertices; + + /*Cleanup*/ + if(segments) xDelete(segments); + xDelete(masklevelset); + xDelete(pelementerror); + +} +/*}}}*/ +void FemModel::InitializeAdaptiveRefinementNeopz(void){/*{{{*/ + + /*Define variables*/ + int my_rank = IssmComm::GetRank(); + int numberofvertices = this->vertices->NumberOfVertices(); + int numberofelements = this->elements->NumberOfElements(); + int numberofsegments = 0; //used on matlab + IssmDouble* x = NULL; + IssmDouble* y = NULL; + IssmDouble* z = NULL; + int* elements = NULL; + int elementswidth = this->GetElementsWidth(); //just tria elements in this version. Itapopo: + int levelmax = 0; + IssmDouble regionlevel1 = 0.; + IssmDouble regionlevelmax = 0.; + + /*Initialize field as NULL for now*/ + this->amr = NULL; + + /*Get vertices coordinates of the coarse mesh (father mesh)*/ + /*elements comes in Matlab indexing*/ + this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elements); + + /*Get amr parameters*/ + this->parameters->FindParam(&levelmax,AmrLevelMaxEnum); + this->parameters->FindParam(®ionlevel1,AmrRegionLevel1Enum); + this->parameters->FindParam(®ionlevelmax,AmrRegionLevelMaxEnum); + + /*Create initial mesh (coarse mesh) in neopz data structure*/ + /*Just CPU #0 should keep AMR object*/ + this->SetRefPatterns(); + if(my_rank==0){ + bool ismisomip = false; + if(ismisomip){//itapopo + TPZFileStream fstr; + std::stringstream ss; + int frictionlaw; + + this->parameters->FindParam(&frictionlaw,FrictionLawEnum); + + ss << levelmax; + if(frictionlaw==1){ + ss << "_viscous/amr.txt"; + }else if(frictionlaw==7){ + ss << "_tsai/amr.txt"; + }else{ + _error_("friction law not supported here."); + } + + std::string AMRfile = "/home/santos/Misomip2/L" + ss.str(); + fstr.OpenRead(AMRfile.c_str()); + + TPZSaveable *sv = TPZSaveable::Restore(fstr,0); + this->amr = dynamic_cast(sv); + } + else{ + this->amr = new AdaptiveMeshRefinement(); + //this->amr->SetLevelMax(levelmax); //Set max level of refinement + //this->amr->SetRegions(regionlevel1,regionlevelmax); + this->amr->CreateInitialMesh(numberofvertices,numberofelements,numberofsegments,elementswidth,x,y,z,elements,NULL); + } + this->amr->SetLevelMax(levelmax); //Set max level of refinement + this->amr->SetRegions(regionlevel1,regionlevelmax); + } + + /*Free the vectors*/ + xDelete(x); + xDelete(y); + xDelete(z); + xDelete(elements); +} +/*}}}*/ +void FemModel::SetRefPatterns(){/*{{{*/ + + /*Initialize the global variable of refinement patterns*/ + gRefDBase.InitializeUniformRefPattern(EOned); + gRefDBase.InitializeUniformRefPattern(ETriangle); + + //gRefDBase.InitializeRefPatterns(); + /*Insert specifics patterns to ISSM core*/ + std::string filepath = REFPATTERNDIR; + std::string filename1 = filepath + "/2D_Triang_Rib_3.rpt"; + std::string filename2 = filepath + "/2D_Triang_Rib_4.rpt"; + std::string filename3 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4.rpt"; + std::string filename4 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4_permuted.rpt"; + std::string filename5 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5.rpt"; + std::string filename6 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5_permuted.rpt"; + std::string filename7 = filepath + "/2D_Triang_Rib_5.rpt"; + + TPZAutoPointer refpat1 = new TPZRefPattern(filename1); + TPZAutoPointer refpat2 = new TPZRefPattern(filename2); + TPZAutoPointer refpat3 = new TPZRefPattern(filename3); + TPZAutoPointer refpat4 = new TPZRefPattern(filename4); + TPZAutoPointer refpat5 = new TPZRefPattern(filename5); + TPZAutoPointer refpat6 = new TPZRefPattern(filename6); + TPZAutoPointer refpat7 = new TPZRefPattern(filename7); + + if(!gRefDBase.FindRefPattern(refpat1)) gRefDBase.InsertRefPattern(refpat1); + if(!gRefDBase.FindRefPattern(refpat2)) gRefDBase.InsertRefPattern(refpat2); + if(!gRefDBase.FindRefPattern(refpat3)) gRefDBase.InsertRefPattern(refpat3); + if(!gRefDBase.FindRefPattern(refpat4)) gRefDBase.InsertRefPattern(refpat4); + if(!gRefDBase.FindRefPattern(refpat5)) gRefDBase.InsertRefPattern(refpat5); + if(!gRefDBase.FindRefPattern(refpat6)) gRefDBase.InsertRefPattern(refpat6); + if(!gRefDBase.FindRefPattern(refpat7)) gRefDBase.InsertRefPattern(refpat7); +} +/*}}}*/ +#endif Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 21801) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 21802) @@ -22,4 +22,7 @@ #ifdef _HAVE_NEOPZ_ #include "./AdaptiveMeshRefinement.h" +#endif +#ifdef _HAVE_BAMG_ +#include "./AmrBamg.h" #endif /*}}}*/ @@ -47,6 +50,11 @@ Vertices *vertices; //one set of vertices + //FIXME: do we want only one class and have virtual functions? or keep 2 classes, at least rename AdaptiveMeshRefinement -> AmrNeopz #ifdef _HAVE_NEOPZ_ AdaptiveMeshRefinement *amr; //adaptive mesh refinement object. It keeps coarse mesh and execute refinement process + #endif + + #ifdef _HAVE_BAMG_ + AmrBamg *amrbamg; //adaptive mesh refinement object. It keeps coarse mesh and execute refinement process #endif @@ -61,7 +69,9 @@ FemModel* copy(); void Echo(); + int GetElementsWidth(){return 3;};//just tria elements in this first version void InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, char* restartfilename, const int solution_type,bool trace,IssmPDouble* X=NULL); void InitFromFids(char* rootpath, FILE* IOMODEL, FILE* toolkitsoptionsfid, int in_solution_type, bool trace, IssmPDouble* X=NULL); void Marshall(char** pmarshalled_data, int* pmarshalled_data_size, int marshall_direction); + void ReMesh(void); void Restart(void); void SetCurrentConfiguration(int configuration_type); @@ -146,13 +156,8 @@ #endif - #ifdef _HAVE_NEOPZ_ - /*Adaptive mesh refinement methods*/ - void InitializeAdaptiveRefinement(void); - void ReMesh(void); + /*AMR*/ void BedrockFromMismipPlus(void); void AdjustBaseThicknessAndMask(void); void GetMesh(Vertices* femmodel_vertices,Elements* femmodel_elements,IssmDouble** px, IssmDouble** py, IssmDouble** pz, int** pelementslist); - int GetElementsWidth(){return 3;};//just tria elements in this first version - void ExecuteRefinement(int &numberofvertices,int &numberofelements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist); void GetGroundediceLevelSet(IssmDouble** pmasklevelset); void CreateVertices(int newnumberofvertices,int newnumberofelements,int elementswidth,int* newelementslist,int* my_vertices,IssmDouble* newx,IssmDouble* newy,IssmDouble* newz,Vertices* vertices); @@ -168,7 +173,16 @@ void SmoothedDeviatoricStressTensor(IssmDouble** ptauxx,IssmDouble** ptauyy,IssmDouble** ptauxy); //nodal values, just for SSA-P1: TauXX, TauYY, TauXY void ZZErrorEstimator(IssmDouble** pelementerror); + + #ifdef _HAVE_BAMG_ + void ReMeshBamg(int* pnewnumberofvertices,int* pnewnumberofelements,IssmDouble** pnewx,IssmDouble** pnewy,IssmDouble** pnewz,int** pnewelementslist); + void InitializeAdaptiveRefinementBamg(void); + #endif + + #ifdef _HAVE_NEOPZ_ + /*Adaptive mesh refinement methods*/ + void ReMeshNeopz(int* pnewnumberofvertices,int* pnewnumberofelements,IssmDouble** pnewx,IssmDouble** pnewy,IssmDouble** pnewz,int** pnewelementslist); + void InitializeAdaptiveRefinementNeopz(void); #endif }; - #endif Index: /issm/trunk-jpl/src/c/classes/IoModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/IoModel.cpp (revision 21801) +++ /issm/trunk-jpl/src/c/classes/IoModel.cpp (revision 21802) @@ -605,4 +605,5 @@ if(strcmp(record_name,"md.materials.type")==0) integer = IoCodeToEnumMaterials(integer); if(strcmp(record_name,"md.timestepping.type")==0) integer = IoCodeToEnumTimestepping(integer); + if(strcmp(record_name,"md.amr.type")==0) integer = IoCodeToEnumAmr(integer); /*Broadcast to other cpus*/ Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 21801) +++ /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 21802) @@ -111,5 +111,4 @@ } } - if(VerboseSolution()) _printf0_(" computing thermal regime\n"); thermal_core(femmodel); } @@ -182,10 +181,16 @@ /*Adaptive mesh refinement*/ - #ifdef _HAVE_NEOPZ_ if(amr_frequency){ if(save_results) femmodel->WriteMeshInResults(); - if(step%amr_frequency==0 && timeReMesh();//Do not refine the last step - } - #endif + if(step%amr_frequency==0 && timeReMesh();//Do not refine the last step + } + //if(save_results){ + // printf("-------------- file: transient_core.cpp line: %i\n",__LINE__); + // femmodel->WriteMeshInResults(); + // OutputResultsx(femmodel); + //} + } } Index: /issm/trunk-jpl/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp (revision 21801) +++ /issm/trunk-jpl/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp (revision 21802) @@ -24,5 +24,5 @@ R2 r; I2 I; - int i,j; + int i,j,k; int it; int i0,i1,i2; @@ -64,4 +64,14 @@ if(y_data[i]ymax) ymax=y_data[i]; + } + + /*Create Single vertex to element connectivity*/ + int* connectivity = xNew(nods_data); + for(i=0;i=0 & k=0 & i0=0 && it=nels_data){ _error_("Triangle number " << it << " not in [0 " << nels_data << "], report bug to developers (interpolation point: " <(reft); + xDelete(connectivity); *pdata_interp=data_interp; return 1; Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 21801) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 21802) @@ -21,4 +21,5 @@ int numoutputs,materialtype,smb_model,basalforcing_model,timestepping_type; char** requestedoutputs = NULL; + char* fieldname = NULL; IssmDouble time; @@ -60,5 +61,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.inversion.type",InversionTypeEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.calving.law",CalvingLawEnum)); - {/*This is specific to ice...*/ + {/*This is specific to ice...*/ parameters->AddObject(iomodel->CopyConstantObject("md.mesh.elementtype",MeshElementtypeEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.steadystate.reltol",SteadystateReltolEnum)); @@ -81,8 +82,4 @@ parameters->AddObject(iomodel->CopyConstantObject("md.materials.rheology_law",MaterialsRheologyLawEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.gia.cross_section_shape",GiaCrossSectionShapeEnum)); - /*amr properties*/ - parameters->AddObject(iomodel->CopyConstantObject("md.amr.level_max",AmrLevelMaxEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.amr.region_level_1",AmrRegionLevel1Enum)); - parameters->AddObject(iomodel->CopyConstantObject("md.amr.region_level_max",AmrRegionLevelMaxEnum)); /*For stress balance only*/ @@ -95,7 +92,38 @@ if(iomodel->domaintype==Domain3DEnum) parameters->AddObject(iomodel->CopyConstantObject("md.mesh.numberoflayers",MeshNumberoflayersEnum)); - } - - + } + + /*amr properties*/ + int amrtype,amr_frequency; + iomodel->FindConstant(&amr_frequency,"md.transient.amr_frequency"); + if(solution_type==TransientSolutionEnum && amr_frequency){ + parameters->AddObject(iomodel->CopyConstantObject("md.amr.type",AmrTypeEnum)); + iomodel->FindConstant(&amrtype,"md.amr.type"); + switch(amrtype){ + #ifdef _HAVE_NEOPZ_ + case AmrNeopzEnum: + parameters->AddObject(iomodel->CopyConstantObject("md.amr.level_max",AmrLevelMaxEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.region_level_1",AmrRegionLevel1Enum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.region_level_max",AmrRegionLevelMaxEnum)); + break; + #endif + + #ifdef _HAVE_BAMG_ + case AmrBamgEnum: + parameters->AddObject(iomodel->CopyConstantObject("md.amr.hmin",AmrHminEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.hmax",AmrHmaxEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.err",AmrErrEnum)); + /*Convert fieldname to enum and put it in params*/ + iomodel->FindConstant(&fieldname,"md.amr.fieldname"); + parameters->AddObject(new IntParam(AmrFieldEnum,StringToEnumx(fieldname))); + xDelete(fieldname); + break; + #endif + + default: + _error_("Adaptive mesh refinement "<AddObject(iomodel->CopyConstantObject("md.basalforcings.model",BasalforcingsEnum)); @@ -255,5 +283,5 @@ ParseToolkitsOptionsx(parameters,toolkitsoptionsfid); - #ifdef _HAVE_ADOLC_ + #ifdef _HAVE_ADOLC_ if(VerboseMProcessor()) _printf0_(" starting autodiff parameters \n"); CreateParametersAutodiff(parameters,iomodel); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21801) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21802) @@ -833,8 +833,15 @@ /*}}}*/ /*Adaptive mesh refinement (AMR){{{*/ + TransientAmrFrequencyEnum, + AmrTypeEnum, + AmrNeopzEnum, AmrLevelMaxEnum, AmrRegionLevel1Enum, AmrRegionLevelMaxEnum, - TransientAmrFrequencyEnum, + AmrBamgEnum, + AmrHminEnum, + AmrHmaxEnum, + AmrFieldEnum, + AmrErrEnum, /*}}}*/ ParametersENDEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21801) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21802) @@ -809,8 +809,15 @@ case EsaRequestedOutputsEnum : return "EsaRequestedOutputs"; case EsaNumRequestedOutputsEnum : return "EsaNumRequestedOutputs"; + case TransientAmrFrequencyEnum : return "TransientAmrFrequency"; + case AmrTypeEnum : return "AmrType"; + case AmrNeopzEnum : return "AmrNeopz"; case AmrLevelMaxEnum : return "AmrLevelMax"; case AmrRegionLevel1Enum : return "AmrRegionLevel1"; case AmrRegionLevelMaxEnum : return "AmrRegionLevelMax"; - case TransientAmrFrequencyEnum : return "TransientAmrFrequency"; + case AmrBamgEnum : return "AmrBamg"; + case AmrHminEnum : return "AmrHmin"; + case AmrHmaxEnum : return "AmrHmax"; + case AmrFieldEnum : return "AmrField"; + case AmrErrEnum : return "AmrErr"; case ParametersENDEnum : return "ParametersEND"; case XYEnum : return "XY"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21801) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21802) @@ -827,8 +827,15 @@ else if (strcmp(name,"EsaRequestedOutputs")==0) return EsaRequestedOutputsEnum; else if (strcmp(name,"EsaNumRequestedOutputs")==0) return EsaNumRequestedOutputsEnum; + else if (strcmp(name,"TransientAmrFrequency")==0) return TransientAmrFrequencyEnum; + else if (strcmp(name,"AmrType")==0) return AmrTypeEnum; + else if (strcmp(name,"AmrNeopz")==0) return AmrNeopzEnum; else if (strcmp(name,"AmrLevelMax")==0) return AmrLevelMaxEnum; else if (strcmp(name,"AmrRegionLevel1")==0) return AmrRegionLevel1Enum; else if (strcmp(name,"AmrRegionLevelMax")==0) return AmrRegionLevelMaxEnum; - else if (strcmp(name,"TransientAmrFrequency")==0) return TransientAmrFrequencyEnum; + else if (strcmp(name,"AmrBamg")==0) return AmrBamgEnum; + else if (strcmp(name,"AmrHmin")==0) return AmrHminEnum; + else if (strcmp(name,"AmrHmax")==0) return AmrHmaxEnum; + else if (strcmp(name,"AmrField")==0) return AmrFieldEnum; + else if (strcmp(name,"AmrErr")==0) return AmrErrEnum; else if (strcmp(name,"ParametersEND")==0) return ParametersENDEnum; else if (strcmp(name,"XY")==0) return XYEnum; @@ -868,5 +875,8 @@ else if (strcmp(name,"Moulin")==0) return MoulinEnum; else if (strcmp(name,"Pengrid")==0) return PengridEnum; - else if (strcmp(name,"Penpair")==0) return PenpairEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"Penpair")==0) return PenpairEnum; else if (strcmp(name,"Profiler")==0) return ProfilerEnum; else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum; @@ -875,8 +885,5 @@ else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum; else if (strcmp(name,"ElementSId")==0) return ElementSIdEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"VectorParam")==0) return VectorParamEnum; + else if (strcmp(name,"VectorParam")==0) return VectorParamEnum; else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum; else if (strcmp(name,"Segment")==0) return SegmentEnum; @@ -991,5 +998,8 @@ else if (strcmp(name,"TotalGroundedBmb")==0) return TotalGroundedBmbEnum; else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum; - else if (strcmp(name,"P0")==0) return P0Enum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"P0")==0) return P0Enum; else if (strcmp(name,"P0Array")==0) return P0ArrayEnum; else if (strcmp(name,"P1")==0) return P1Enum; @@ -998,8 +1008,5 @@ else if (strcmp(name,"P1bubblecondensed")==0) return P1bubblecondensedEnum; else if (strcmp(name,"P2")==0) return P2Enum; - else stage=9; - } - if(stage==9){ - if (strcmp(name,"P2bubble")==0) return P2bubbleEnum; + else if (strcmp(name,"P2bubble")==0) return P2bubbleEnum; else if (strcmp(name,"P2bubblecondensed")==0) return P2bubblecondensedEnum; else if (strcmp(name,"P2xP1")==0) return P2xP1Enum; Index: /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 21801) +++ /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 21802) @@ -128,4 +128,11 @@ } }/*}}}*/ +int IoCodeToEnumAmr(int enum_in){/*{{{*/ + switch(enum_in){ + case 1: return AmrBamgEnum; + case 2: return AmrNeopzEnum; + default: _error_("Marshalled AMR code \""<',0,'NaN',1); + md = checkfield(md,'fieldname','amr.hmin','numel',[1],'>',0,'<',self.hmax,'NaN',1); + %md = checkfield(md,'fieldname','amr.fieldname','string',[1]); + end % }}} + function disp(self) % {{{ + disp(sprintf(' amrbamg parameters:')); + + fielddisplay(self,'hmin',['minimum element length']); + fielddisplay(self,'hmax',['maximum element length']); + fielddisplay(self,'fieldname',['name of input that will be used to compute the metric (should be an input of FemModel)']); + + end % }}} + function marshall(self,prefix,md,fid) % {{{ + + WriteData(fid,prefix,'name','md.amr.type','data',1,'format','Integer'); + WriteData(fid,prefix,'object',self,'class','amr','fieldname','hmin','format','Double'); + WriteData(fid,prefix,'object',self,'class','amr','fieldname','hmax','format','Double'); + WriteData(fid,prefix,'object',self,'class','amr','fieldname','fieldname','format','String'); + WriteData(fid,prefix,'object',self,'class','err','fieldname','err','format','Double'); + + end % }}} + end +end Index: /issm/trunk-jpl/src/m/classes/trans.js =================================================================== --- /issm/trunk-jpl/src/m/classes/trans.js (revision 21801) +++ /issm/trunk-jpl/src/m/classes/trans.js (revision 21802) @@ -9,16 +9,16 @@ //full analysis: Stressbalance, Masstransport and Thermal but no groundingline migration for now - this.issmb = 1; - this.ismasstransport = 1; - this.isstressbalance = 1; - this.isthermal = 1; - this.isgroundingline = 0; - this.isgia = 0; + this.issmb = 1; + this.ismasstransport = 1; + this.isstressbalance = 1; + this.isthermal = 1; + this.isgroundingline = 0; + this.isgia = 0; this.isdamageevolution = 0; - this.ismovingfront = 0; - this.ishydrology = 0; - this.isslr = 0; - this.iscoupler = 0; - this.amr_frequency = 1; + this.ismovingfront = 0; + this.ishydrology = 0; + this.isslr = 0; + this.iscoupler = 0; + this.amr_frequency = 0; //default output Index: /issm/trunk-jpl/src/m/classes/transient.m =================================================================== --- /issm/trunk-jpl/src/m/classes/transient.m (revision 21801) +++ /issm/trunk-jpl/src/m/classes/transient.m (revision 21802) @@ -68,5 +68,5 @@ self.isoceancoupling = 0; self.iscoupler = 0; - self.amr_frequency = 1; + self.amr_frequency = 0; %default output Index: /issm/trunk-jpl/src/m/classes/transient.py =================================================================== --- /issm/trunk-jpl/src/m/classes/transient.py (revision 21801) +++ /issm/trunk-jpl/src/m/classes/transient.py (revision 21802) @@ -84,18 +84,18 @@ #full analysis: Stressbalance, Masstransport and Thermal but no groundingline migration for now - self.issmb = True - self.ismasstransport = True - self.isstressbalance = True - self.isthermal = True - self.isgroundingline = False - self.isgia = False - self.isesa = False + self.issmb = True + self.ismasstransport = True + self.isstressbalance = True + self.isthermal = True + self.isgroundingline = False + self.isgia = False + self.isesa = False self.isdamageevolution = False - self.ismovingfront = False - self.ishydrology = False - self.isslr = False - self.isoceancoupling = False - self.iscoupler = False - self.amr_frequency = 1 + self.ismovingfront = False + self.ishydrology = False + self.isslr = False + self.isoceancoupling = False + self.iscoupler = False + self.amr_frequency = 0 #default output