![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
source:
issm/oecreview/Archive/22819-23185/ISSM-23045-23046.diff
| Last change on this file was 23186, checked in by , 7 years ago | |
|---|---|
| File size: 62.9 KB | |
-
../trunk-jpl/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp
1 1 /*!\file: InterpFromGridToMeshx.cpp 2 2 * \brief "c" core code for interpolating values from a structured grid. 3 */ 3 */ 4 4 5 5 #ifdef HAVE_CONFIG_H 6 6 #include <config.h> … … 24 24 double* y=NULL; 25 25 int i; 26 26 27 // TEST28 _printf_("\r N: "<<N<<" M:"<<M<<"\n");29 30 27 /*Some checks on arguments: */ 31 28 if ((M<2) || (N<2) || (nods<=0)){ 32 29 _error_("nothing to be done according to the dimensions of input matrices and vectors."); … … 150 147 * | | | 151 148 * | | | 152 149 * y1 x---------+-----x Q21 153 * x1 x2 150 * x1 x2 154 151 * 155 152 */ 156 153 x1=x[n]; x2=x[n+1]; … … 224 221 double triangleinterp(double x1,double x2,double y1,double y2,double Q11,double Q12,double Q21,double Q22,double x,double y){ 225 222 /*split the rectangle in 2 triangle and 226 223 * use Lagrange P1 interpolation 227 * 224 * 228 225 * +3----+2,3' Q12----Q22 229 226 * | /| | /| 230 227 * | / | | / | … … 281 278 _assert_(x2>x1 && y2>y1); 282 279 _assert_(x<=x2 && x>=x1 && y<=y2 && y>=y1); 283 280 284 return 281 return 285 282 +Q11*(x2-x)*(y2-y)/((x2-x1)*(y2-y1)) 286 283 +Q21*(x-x1)*(y2-y)/((x2-x1)*(y2-y1)) 287 284 +Q12*(x2-x)*(y-y1)/((x2-x1)*(y2-y1)) … … 300 297 * | | | 301 298 * | | | 302 299 * y1 x--------x---------x Q21 303 * x1 xm x2 300 * x1 xm x2 304 301 * 305 302 */ 306 303 /*Checks*/ -
../trunk-jpl/src/wrappers/InterpFromGridToMesh/InterpFromGridToMesh.js
5 5 * 6 6 * Usage: 7 7 * var data_mesh=InterpFromGridToMesh(xIn,yIn,dataIn,xMeshIn,yMeshIn,defaultValue,interpType);\ 8 * 8 * 9 9 * xIn,yIn : coordinates of matrix data. (x and y must be in increasing order) 10 10 * dataIn : matrix holding the data to be interpolated onto the mesh 11 11 * xMeshIn,yMeshIn : coordinates of the points onto which we interpolate … … 50 50 var y = {}; 51 51 var yMesh = {}; 52 52 //}}} 53 54 53 54 55 55 /* 56 56 Dynamic allocations 57 57 */ 58 58 //{{{ 59 59 60 60 /* 61 61 Input 62 62 */ … … 67 67 dxHeap = new Uint8Array(Module.HEAPU8.buffer, dxPtr, nx); 68 68 dxHeap.set(new Uint8Array(dx.buffer)); 69 69 x = dxHeap.byteOffset; 70 70 71 71 dy = new Float64Array(yIn); 72 72 ny = dy.length * dy.BYTES_PER_ELEMENT; 73 73 dyPtr = Module._malloc(ny); … … 74 74 dyHeap = new Uint8Array(Module.HEAPU8.buffer, dyPtr, ny); 75 75 dyHeap.set(new Uint8Array(dy.buffer)); 76 76 y = dyHeap.byteOffset; 77 77 78 78 ddata = new Float64Array(dataIn); 79 79 ndata = ddata.length * ddata.BYTES_PER_ELEMENT; 80 80 ddataPtr = Module._malloc(ndata); … … 81 81 ddataHeap = new Uint8Array(Module.HEAPU8.buffer, ddataPtr, ndata); 82 82 ddataHeap.set(new Uint8Array(ddata.buffer)); 83 83 data = ddataHeap.byteOffset; 84 84 85 85 dxMesh = new Float64Array(xMeshIn); 86 86 nxMesh = dxMesh.length * dxMesh.BYTES_PER_ELEMENT; 87 87 dxMeshPtr = Module._malloc(nxMesh); … … 88 88 dxMeshHeap = new Uint8Array(Module.HEAPU8.buffer, dxMeshPtr, nxMesh); 89 89 dxMeshHeap.set(new Uint8Array(dxMesh.buffer)); 90 90 xMesh = dxMeshHeap.byteOffset; 91 91 92 92 dyMesh = new Float64Array(yMeshIn); 93 93 nyMesh = dyMesh.length * dyMesh.BYTES_PER_ELEMENT; 94 94 dyMeshPtr = Module._malloc(nyMesh); … … 95 95 dyMeshHeap = new Uint8Array(Module.HEAPU8.buffer, dyMeshPtr, nyMesh); 96 96 dyMeshHeap.set(new Uint8Array(dyMesh.buffer)); 97 97 yMesh = dyMeshHeap.byteOffset; 98 98 99 99 nods = xMeshIn.length; 100 100 meshNumRows = xMeshIn.length; 101 console.log('InterpFromGridToMesh.js: meshNumRows => ' + meshNumRows); 102 103 101 102 104 103 /* 105 104 Retrieve interpolation type 106 105 */ … … 111 110 interpType = 'bilinear'; 112 111 } 113 112 //}}} 114 113 115 114 /* 116 115 Output 117 116 */ 118 117 pdataMesh = Module._malloc(4); 119 118 //}}} 120 119 121 120 //}}} 122 123 121 122 124 123 /* 125 124 Declare InterpFromGridToMesh module 126 125 */ 127 126 //{{{ 128 127 InterpFromGridToMeshModule = Module.cwrap( 129 'InterpFromGridToMeshModule', 130 'number', 128 'InterpFromGridToMeshModule', 129 'number', 131 130 [ 132 131 'number', // output : pdataMesh 133 132 'number', // input : x … … 144 143 ] 145 144 ); 146 145 //}}} 147 148 146 147 149 148 /* 150 149 Call InterpFromGridToMesh module 151 150 */ 152 151 //{{{ 153 152 InterpFromGridToMeshModule( 154 pdataMesh, 155 x, 156 y, 157 data, 158 xMesh, 159 yMesh, 160 defaultValue, 153 pdataMesh, 154 x, 155 y, 156 data, 157 xMesh, 158 yMesh, 159 defaultValue, 161 160 nods, 162 161 dataNumRowsIn, 163 162 dataNumColsIn, … … 165 164 interpType 166 165 ); 167 166 //}}} 168 169 167 168 170 169 /* 171 170 Dynamic copying from heap 172 171 */ … … 174 173 dataMeshPtr = Module.getValue(pdataMesh, 'i32'); 175 174 dataMesh = Module.HEAPF64.slice(dataMeshPtr / 8, dataMeshPtr / 8 + nods); 176 175 //}}} 177 178 176 177 179 178 /* 180 179 Free resources 181 180 */ … … 182 181 //{{{ 183 182 Module._free(pdataMesh); 184 183 //}}} 185 186 184 185 187 186 return dataMesh; 188 187 } -
../trunk-jpl/jenkins/macosx_pine-island
4 4 #-------------------------------# 5 5 6 6 #MATLAB path 7 MATLAB_PATH="/Applications/MATLAB_R201 7b.app/"7 MATLAB_PATH="/Applications/MATLAB_R2015b.app/" 8 8 9 #ISSM CONFIGURATION 9 #ISSM CONFIGURATION 10 10 ISSM_CONFIG='--prefix=$ISSM_DIR \ 11 11 --with-matlab-dir=$MATLAB_PATH \ 12 12 --with-triangle-dir=$ISSM_DIR/externalpackages/triangle/install \ … … 53 53 #as follows: runme($MATLAB_NROPTIONS). The options must be understandable 54 54 #by Matlab and runme.m 55 55 #ex: "'id',[101 102 103]" 56 ## FS 56 ## FS 57 57 PYTHON_NROPTIONS="" 58 58 MATLAB_NROPTIONS="'exclude',[701,702,703,435,IdFromString('Dakota')]" -
../trunk-jpl/jenkins/macosx_pine-island_static
4 4 #-------------------------------# 5 5 6 6 #MATLAB path 7 MATLAB_PATH="/Applications/MATLAB_R201 7b.app/"7 MATLAB_PATH="/Applications/MATLAB_R2015b.app/" 8 8 9 #ISSM CONFIGURATION 9 #ISSM CONFIGURATION 10 10 ISSM_CONFIG='--prefix=$ISSM_DIR \ 11 11 --disable-static \ 12 12 --enable-standalone-executables \ … … 22 22 --with-metis-dir=$ISSM_DIR/externalpackages/petsc/install \ 23 23 --with-m1qn3-dir=$ISSM_DIR/externalpackages/m1qn3/install \ 24 24 --with-math77-dir=$ISSM_DIR/externalpackages/math77/install \ 25 --with-fortran-lib="/usr/local/gfortran/lib/libgfortran.a /usr/local/gfortran/lib/libquadmath.a /usr/local/gfortran/lib/gcc/x86_64-apple-darwin1 5/6.1.0/libgcc.a" \25 --with-fortran-lib="/usr/local/gfortran/lib/libgfortran.a /usr/local/gfortran/lib/libquadmath.a /usr/local/gfortran/lib/gcc/x86_64-apple-darwin14/5.2.0/libgcc.a" \ 26 26 --with-numthreads=4' 27 27 28 28 #PYTHON and MATLAB testing … … 60 60 #as follows: runme($MATLAB_NROPTIONS). The options must be understandable 61 61 #by Matlab and runme.m 62 62 #ex: "'id',[101 102 103]" 63 ## bamg mesh FS 63 ## bamg mesh FS 64 64 #PYTHON_NROPTIONS="" 65 65 #MATLAB_NROPTIONS="'exclude',[119,243,514,701,702,703,435,IdFromString('Dakota')]" -
../trunk-jpl/src/c/classes/FemModel.cpp
1463 1463 1464 1464 /*Figure out maximum across the cluster: */ 1465 1465 ISSM_MPI_Reduce(&maxvy,&node_maxvy,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); 1466 ISSM_MPI_Bcast(&node_maxvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 1466 ISSM_MPI_Bcast(&node_maxvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 1467 1467 maxvy=node_maxvy; 1468 1468 1469 1469 /*Assign output pointers:*/ … … 1487 1487 1488 1488 /*Figure out maximum across the cluster: */ 1489 1489 ISSM_MPI_Reduce(&maxvz,&node_maxvz,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); 1490 ISSM_MPI_Bcast(&node_maxvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 1490 ISSM_MPI_Bcast(&node_maxvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 1491 1491 maxvz=node_maxvz; 1492 1492 1493 1493 /*Assign output pointers:*/ … … 1511 1511 1512 1512 /*Figure out minimum across the cluster: */ 1513 1513 ISSM_MPI_Reduce(&minvel,&node_minvel,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); 1514 ISSM_MPI_Bcast(&node_minvel,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 1514 ISSM_MPI_Bcast(&node_minvel,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 1515 1515 minvel=node_minvel; 1516 1516 1517 1517 /*Assign output pointers:*/ … … 1535 1535 1536 1536 /*Figure out minimum across the cluster: */ 1537 1537 ISSM_MPI_Reduce(&minvx,&node_minvx,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); 1538 ISSM_MPI_Bcast(&node_minvx,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 1538 ISSM_MPI_Bcast(&node_minvx,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 1539 1539 minvx=node_minvx; 1540 1540 1541 1541 /*Assign output pointers:*/ … … 1559 1559 1560 1560 /*Figure out minimum across the cluster: */ 1561 1561 ISSM_MPI_Reduce(&minvy,&node_minvy,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); 1562 ISSM_MPI_Bcast(&node_minvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 1562 ISSM_MPI_Bcast(&node_minvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 1563 1563 minvy=node_minvy; 1564 1564 1565 1565 /*Assign output pointers:*/ … … 1583 1583 1584 1584 /*Figure out minimum across the cluster: */ 1585 1585 ISSM_MPI_Reduce(&minvz,&node_minvz,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); 1586 ISSM_MPI_Bcast(&node_minvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 1586 ISSM_MPI_Bcast(&node_minvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 1587 1587 minvz=node_minvz; 1588 1588 1589 1589 /*Assign output pointers:*/ … … 1628 1628 weights_input->GetInputValue(&weight,gauss,OmegaAbsGradientEnum); 1629 1629 omega_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss); 1630 1630 1631 /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ 1631 /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ 1632 1632 //J+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight; 1633 1633 J+=weight*1/2*pow(dp[0]*dp[0]+dp[1]*dp[1],2)*Jdet*gauss->weight; 1634 1634 } … … 1964 1964 else{ 1965 1965 IssmDouble* values = xNewZeroInit<IssmDouble>(nlines*ncols); 1966 1966 IssmDouble* allvalues = xNew<IssmDouble>(nlines*ncols); 1967 1967 1968 1968 /*Fill-in matrix*/ 1969 1969 for(int j=0;j<elements->Size();j++){ 1970 1970 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(j)); … … 1973 1973 /*Gather from all cpus*/ 1974 1974 ISSM_MPI_Allreduce((void*)values,(void*)allvalues,ncols*nlines,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); 1975 1975 xDelete<IssmDouble>(values); 1976 1976 1977 1977 if(save_results)results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,output_enum,allvalues,nlines,ncols,step,time)); 1978 1978 xDelete<IssmDouble>(allvalues); 1979 1979 } … … 2074 2074 case MaterialsRheologyBbarEnum: this->ElementResponsex(responses,MaterialsRheologyBbarEnum); break; 2075 2075 case VelEnum: this->ElementResponsex(responses,VelEnum); break; 2076 2076 case FrictionCoefficientEnum: NodalValuex(responses, FrictionCoefficientEnum,elements,nodes, vertices, loads, materials, parameters); break; 2077 default: 2077 default: 2078 2078 if(response_descriptor_enum>=Outputdefinition1Enum && response_descriptor_enum <=Outputdefinition100Enum){ 2079 2079 IssmDouble double_result = OutputDefinitionsResponsex(this,response_descriptor_enum); 2080 2080 *responses=double_result; 2081 2081 } 2082 else _error_("response descriptor \"" << EnumToStringx(response_descriptor_enum) << "\" not supported yet!"); 2083 break; 2082 else _error_("response descriptor \"" << EnumToStringx(response_descriptor_enum) << "\" not supported yet!"); 2083 break; 2084 2084 } 2085 2085 2086 2086 } … … 2211 2211 weights_input->GetInputValue(&weight,gauss,ThicknessAbsGradientEnum); 2212 2212 thickness_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss); 2213 2213 2214 /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ 2214 /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ 2215 2215 J+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight; 2216 2216 } 2217 2217 … … 2402 2402 Analysis* analysis= EnumToAnalysis(analysis_type); 2403 2403 analysis->UpdateConstraints(this); 2404 2404 delete analysis; 2405 2405 2406 2406 /*Second, constraints might be time dependent: */ 2407 SpcNodesx(nodes,constraints,parameters,analysis_type); 2407 SpcNodesx(nodes,constraints,parameters,analysis_type); 2408 2408 2409 2409 /*Now, update degrees of freedoms: */ 2410 2410 NodesDofx(nodes,parameters,analysis_type); … … 2415 2415 2416 2416 IssmDouble *surface = NULL; 2417 2417 IssmDouble *bed = NULL; 2418 2418 2419 2419 if(VerboseSolution()) _printf0_(" updating vertices positions\n"); 2420 2420 2421 2421 /*get vertex vectors for bed and thickness: */ … … 2454 2454 /*}}}*/ 2455 2455 2456 2456 /*AMR*/ 2457 #if !defined(_HAVE_ADOLC_) 2457 #if !defined(_HAVE_ADOLC_) 2458 2458 void FemModel::ReMesh(void){/*{{{*/ 2459 2459 2460 2460 /*Intermediaries*/ … … 2464 2464 int *newelementslist = NULL; 2465 2465 int newnumberofvertices = -1; 2466 2466 int newnumberofelements = -1; 2467 bool* my_elements = NULL; 2467 bool* my_elements = NULL; 2468 2468 int* my_vertices = NULL; 2469 2469 int elementswidth = this->GetElementsWidth();//just tria elements in this version 2470 2470 int amrtype,basalforcing_model; 2471 2471 bool isgroundingline; 2472 2472 2473 /*Branch to specific amr depending on requested method*/ 2473 /*Branch to specific amr depending on requested method*/ 2474 2474 this->parameters->FindParam(&amrtype,AmrTypeEnum); 2475 2475 switch(amrtype){ 2476 2476 #if defined(_HAVE_NEOPZ_) && !defined(_HAVE_ADOLC_) … … 2483 2483 2484 2484 default: _error_("not implemented yet"); 2485 2485 } 2486 2486 2487 2487 /*Partitioning the new mesh. Maybe ElementsAndVerticesPartitioning.cpp could be modified to set this without iomodel.*/ 2488 2488 this->ElementsAndVerticesPartitioning(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,&my_elements,&my_vertices); 2489 2489 … … 2514 2514 int analysis_enum = this->analysis_type_list[i]; 2515 2515 2516 2516 /*As the domain is 2D, it is not necessary to create nodes for this analysis*/ 2517 if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue; 2517 if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue; 2518 2518 2519 2519 this->CreateNodes(newnumberofvertices,my_vertices,nodecounter,analysis_enum,new_nodes); 2520 2520 this->CreateConstraints(new_vertices,nodecounter,constraintcounter,analysis_enum,new_constraints); … … 2544 2544 analysis_type=this->analysis_type_list[i]; 2545 2545 //SetCurrentConfiguration(analysis_type); 2546 2546 2547 this->analysis_counter=i; 2547 this->analysis_counter=i; 2548 2548 /*Now, plug analysis_counter and analysis_type inside the parameters: */ 2549 2549 this->parameters->SetParam(this->analysis_counter,AnalysisCounterEnum); 2550 2550 this->parameters->SetParam(analysis_type,AnalysisTypeEnum); … … 2561 2561 } 2562 2562 2563 2563 ConfigureObjectsx(new_elements,this->loads,new_nodes,new_vertices,new_materials,this->parameters); 2564 if(i==0){ 2564 if(i==0){ 2565 2565 VerticesDofx(new_vertices,this->parameters); //only call once, we only have one set of vertices 2566 2566 } 2567 2567 SpcNodesx(new_nodes,new_constraints,this->parameters,analysis_type); … … 2605 2605 2606 2606 /*Insert bedrock from mismip+ setup*/ 2607 2607 /*This was used to Misomip project/simulations*/ 2608 2608 2609 2609 if(VerboseSolution())_printf0_(" call Mismip bedrock adjust module\n"); 2610 2610 2611 2611 IssmDouble x,y,bx,by; … … 2622 2622 bx = -150.-728.8*pow(x/300000.,2)+343.91*pow(x/300000.,4)-50.57*pow(x/300000.,6); 2623 2623 by = 500./(1.+exp((-2./4000.)*(y-80000./2.-24000.)))+500./(1.+exp((2./4000.)*(y-80000./2.+24000.))); 2624 2624 r[i] = max(bx+by,-720.); 2625 } 2625 } 2626 2626 /*insert new bedrock*/ 2627 2627 element->AddInput(BedEnum,&r[0],P1Enum); 2628 2628 /*Cleanup*/ … … 2635 2635 void FemModel::AdjustBaseThicknessAndMask(void){/*{{{*/ 2636 2636 2637 2637 if(VerboseSolution())_printf0_(" call adjust base and thickness module\n"); 2638 2638 2639 2639 int numvertices = this->GetElementsWidth(); 2640 2640 IssmDouble rho_water,rho_ice,density,base_float; 2641 2641 IssmDouble* phi = xNew<IssmDouble>(numvertices); … … 2647 2647 2648 2648 for(int el=0;el<this->elements->Size();el++){ 2649 2649 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(el)); 2650 2650 2651 2651 element->GetInputListOnVertices(&s[0],SurfaceEnum); 2652 2652 element->GetInputListOnVertices(&r[0],BedEnum); 2653 2653 element->GetInputListOnVertices(&sl[0],SealevelEnum); … … 2659 2659 /*calculate base floatation (which supports given surface*/ 2660 2660 base_float = rho_ice*s[i]/(rho_ice-rho_water); 2661 2661 if(r[i]>base_float){ 2662 b[i] = r[i]; 2663 } 2662 b[i] = r[i]; 2663 } 2664 2664 else { 2665 b[i] = base_float; 2666 } 2665 b[i] = base_float; 2666 } 2667 2667 2668 2668 if(fabs(sl[i])>0) _error_("Sea level value "<<sl[i]<<" not supported!"); 2669 2669 /*update thickness and mask grounded ice level set*/ 2670 2670 h[i] = s[i]-b[i]; 2671 phi[i] = h[i]+r[i]/density; 2671 phi[i] = h[i]+r[i]/density; 2672 2672 } 2673 2673 2674 2674 /*Update inputs*/ … … 2676 2676 element->AddInput(ThicknessEnum,&h[0],P1Enum); 2677 2677 element->AddInput(BaseEnum,&b[0],P1Enum); 2678 2678 } 2679 2679 2680 2680 /*Delete*/ 2681 2681 xDelete<IssmDouble>(phi); 2682 2682 xDelete<IssmDouble>(h); … … 2704 2704 int* P1input_interp = NULL; 2705 2705 Vector<IssmDouble>* input_interpolations = NULL; 2706 2706 IssmDouble* input_interpolations_serial = NULL; 2707 int* pos = NULL; 2707 int* pos = NULL; 2708 2708 IssmDouble value = 0; 2709 2709 2710 2710 /*Figure out how many inputs we have and their respective interpolation*/ … … 2724 2724 P0input_enums = xNew<int>(numP0inputs); 2725 2725 P0input_interp = xNew<int>(numP0inputs); 2726 2726 P1input_enums = xNew<int>(numP1inputs); 2727 P1input_interp = xNew<int>(numP1inputs); 2727 P1input_interp = xNew<int>(numP1inputs); 2728 2728 } 2729 2729 numP0inputs = 0; 2730 2730 numP1inputs = 0; … … 2756 2756 } 2757 2757 } 2758 2758 } 2759 2760 /*Get P0 and P1 inputs over the elements*/ 2759 2760 /*Get P0 and P1 inputs over the elements*/ 2761 2761 pos = xNew<int>(elementswidth); 2762 2762 vP0inputs= new Vector<IssmDouble>(numberofelements*numP0inputs); 2763 2763 vP1inputs= new Vector<IssmDouble>(numberofvertices*numP1inputs); 2764 2764 for(int i=0;i<this->elements->Size();i++){ 2765 2765 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i)); 2766 2766 2767 2767 /*Get P0 inputs*/ 2768 2768 for(int j=0;j<numP0inputs;j++){ 2769 TriaInput* input=xDynamicCast<TriaInput*>(element->GetInput(P0input_enums[j])); 2769 TriaInput* input=xDynamicCast<TriaInput*>(element->GetInput(P0input_enums[j])); 2770 2770 input->GetInputAverage(&value); 2771 2771 pos[0]=element->Sid()*numP0inputs+j; 2772 /*Insert input in the vector*/ 2772 /*Insert input in the vector*/ 2773 2773 vP0inputs->SetValues(1,pos,&value,INS_VAL); 2774 2774 } 2775 2775 2776 2776 /*Get P1 inputs*/ 2777 2777 for(int j=0;j<numP1inputs;j++){ 2778 2778 TriaInput* input=xDynamicCast<TriaInput*>(element->GetInput(P1input_enums[j])); … … 2779 2779 pos[0]=element->vertices[0]->Sid()*numP1inputs+j; 2780 2780 pos[1]=element->vertices[1]->Sid()*numP1inputs+j; 2781 2781 pos[2]=element->vertices[2]->Sid()*numP1inputs+j; 2782 /*Insert input in the vector*/ 2783 vP1inputs->SetValues(elementswidth,pos,input->values,INS_VAL); 2782 /*Insert input in the vector*/ 2783 vP1inputs->SetValues(elementswidth,pos,input->values,INS_VAL); 2784 2784 } 2785 2785 } 2786 2786 … … 2789 2789 vP1inputs->Assemble(); 2790 2790 P0inputs=vP0inputs->ToMPISerial(); 2791 2791 P1inputs=vP1inputs->ToMPISerial(); 2792 2792 2793 2793 /*Assign pointers*/ 2794 *pnumP0inputs = numP0inputs; 2795 *pP0inputs = P0inputs; 2794 *pnumP0inputs = numP0inputs; 2795 *pP0inputs = P0inputs; 2796 2796 *pP0input_enums = P0input_enums; 2797 2797 *pP0input_interp = P0input_interp; 2798 *pnumP1inputs = numP1inputs; 2799 *pP1inputs = P1inputs; 2798 *pnumP1inputs = numP1inputs; 2799 *pP1inputs = P1inputs; 2800 2800 *pP1input_enums = P1input_enums; 2801 *pP1input_interp = P1input_interp; 2801 *pP1input_interp = P1input_interp; 2802 2802 2803 2803 /*Cleanup*/ 2804 2804 delete input_interpolations; … … 2809 2809 } 2810 2810 /*}}}*/ 2811 2811 void FemModel::InterpolateInputs(Vertices* newfemmodel_vertices,Elements* newfemmodel_elements){/*{{{*/ 2812 2812 2813 2813 int numberofelements = this->elements->NumberOfElements(); //global, entire old mesh 2814 2814 int newnumberofelements = newfemmodel_elements->Size(); //just on the new partition 2815 2815 int numberofvertices = this->vertices->NumberOfVertices(); //global, entire old mesh … … 2825 2825 IssmDouble* newP1inputs = NULL; //just on the new partition 2826 2826 int* P1input_enums = NULL; 2827 2827 int* P1input_interp = NULL; 2828 IssmDouble* values = NULL; 2828 IssmDouble* values = NULL; 2829 2829 IssmDouble* vector = NULL; 2830 2830 IssmDouble* x = NULL;//global, entire old mesh 2831 2831 IssmDouble* y = NULL;//global, entire old mesh … … 2837 2837 IssmDouble* newxc = NULL;//just on the new partition 2838 2838 IssmDouble* newyc = NULL;//just on the new partition 2839 2839 int* newelementslist = NULL;//just on the new partition 2840 int* sidtoindex = NULL;//global vertices sid to partition index 2840 int* sidtoindex = NULL;//global vertices sid to partition index 2841 2841 2842 2842 /*Get old P0 and P1 inputs (entire mesh)*/ 2843 2843 this->GetInputs(&numP0inputs,&P0inputs,&P0input_enums,&P0input_interp,&numP1inputs,&P1inputs,&P1input_enums,&P1input_interp); … … 2870 2870 newx,newy,newnumberofvertices,NULL); 2871 2871 2872 2872 /*Insert P0 and P1 inputs into the new elements (just on the new partition)*/ 2873 values=xNew<IssmDouble>(elementswidth); 2873 values=xNew<IssmDouble>(elementswidth); 2874 2874 for(int i=0;i<newfemmodel_elements->Size();i++){//just on the new partition 2875 2875 Element* element=xDynamicCast<Element*>(newfemmodel_elements->GetObjectByOffset(i)); 2876 2876 /*newP0inputs is just on the new partition*/ 2877 2877 for(int j=0;j<numP0inputs;j++){ 2878 switch(P0input_interp[j]){ 2878 switch(P0input_interp[j]){ 2879 2879 case P0Enum: 2880 2880 case DoubleInputEnum: 2881 2881 element->AddInput(new DoubleInput(P0input_enums[j],newP0inputs[i*numP0inputs+j])); 2882 2882 break; 2883 case IntInputEnum: 2883 case IntInputEnum: 2884 2884 element->AddInput(new IntInput(P0input_enums[j],reCast<int>(newP0inputs[i*numP0inputs+j]))); 2885 2885 break; 2886 2886 case BoolInputEnum: … … 2898 2898 element->inputs->AddInput(new TriaInput(P1input_enums[j],values,P1Enum)); 2899 2899 } 2900 2900 } 2901 2901 2902 2902 /*Cleanup*/ 2903 2903 xDelete<IssmDouble>(P0inputs); 2904 2904 xDelete<IssmDouble>(newP0inputs); … … 2937 2937 int* elementslist = NULL; 2938 2938 2939 2939 if(!this->elements || !this->vertices || !this->results || !this->parameters) return; 2940 2940 2941 2941 parameters->FindParam(&step,StepEnum); 2942 2942 parameters->FindParam(&time,TimeEnum); 2943 2943 numberofelements=this->elements->NumberOfElements(); … … 2955 2955 2956 2956 this->results->AddResult(new GenericExternalResult<IssmDouble*>(this->results->Size()+1,MeshYEnum, 2957 2957 y,numberofvertices,1,step,time)); 2958 2958 2959 2959 /*Cleanup*/ 2960 2960 xDelete<IssmDouble>(x); 2961 2961 xDelete<IssmDouble>(y); … … 3016 3016 3017 3017 /*Assemble*/ 3018 3018 vmasklevelset->Assemble(); 3019 3019 3020 3020 /*Serialize and set output*/ 3021 3021 (*pmasklevelset)=vmasklevelset->ToMPISerial(); 3022 3022 … … 3030 3030 void FemModel::CreateVertices(int newnumberofvertices,int newnumberofelements,int elementswidth,int* newelementslist,int* my_vertices,IssmDouble* newx,IssmDouble* newy,IssmDouble* newz,Vertices* vertices){/*{{{*/ 3031 3031 3032 3032 /*newelementslist is in Matlab indexing*/ 3033 3033 3034 3034 /*Creating connectivity table*/ 3035 3035 int* connectivity=NULL; 3036 3036 connectivity=xNewZeroInit<int>(newnumberofvertices); … … 3041 3041 _assert_(vertexid>0 && vertexid-1<newnumberofvertices);//Matlab indexing 3042 3042 connectivity[vertexid-1]+=1;//Matlab to C indexing 3043 3043 } 3044 } 3044 } 3045 3045 3046 3046 /*Create vertex and insert in vertices*/ 3047 3047 for(int i=0;i<newnumberofvertices;i++){ 3048 3048 if(my_vertices[i]){ 3049 Vertex *newvertex=new Vertex(); 3049 Vertex *newvertex=new Vertex(); 3050 3050 newvertex->id=i+1; 3051 3051 newvertex->sid=i; 3052 3052 newvertex->pid=UNDEF; … … 3057 3057 newvertex->sigma=0.; 3058 3058 newvertex->connectivity=connectivity[i]; 3059 3059 newvertex->clone=false;//itapopo check this 3060 vertices->AddObject(newvertex); 3061 } 3060 vertices->AddObject(newvertex); 3061 } 3062 3062 } 3063 3063 3064 3064 xDelete<int>(connectivity); … … 3085 3085 for(int j=0;j<nummodels;j++) newtria->element_type_list[j]=0; 3086 3086 } 3087 3087 else newtria->element_type_list=NULL; 3088 3088 3089 3089 /*Element hook*/ 3090 3090 int matpar_id=newnumberofelements+1; //retrieve material parameter id (last pointer in femodel->materials) 3091 3091 int material_id=i+1; // retrieve material_id = i+1; 3092 3092 /*retrieve vertices ids*/ 3093 3093 int* vertex_ids=xNew<int>(elementswidth); 3094 for(int j=0;j<elementswidth;j++) vertex_ids[j]=reCast<int>(newelementslist[elementswidth*i+j]);//this Hook wants Matlab indexing 3094 for(int j=0;j<elementswidth;j++) vertex_ids[j]=reCast<int>(newelementslist[elementswidth*i+j]);//this Hook wants Matlab indexing 3095 3095 /*Setting the hooks*/ 3096 3096 newtria->numanalyses =this->nummodels; 3097 3097 newtria->hnodes =new Hook*[this->nummodels]; … … 3103 3103 for(int j=0;j<this->nummodels;j++) newtria->hnodes[j]=NULL; 3104 3104 /*Clean up*/ 3105 3105 xDelete<int>(vertex_ids); 3106 elements->AddObject(newtria); 3107 } 3106 elements->AddObject(newtria); 3107 } 3108 3108 } 3109 3109 3110 3110 } … … 3114 3114 /*Just Matice in this version*/ 3115 3115 for(int i=0;i<newnumberofelements;i++){ 3116 3116 if(my_elements[i]){ 3117 materials->AddObject(new Matice(i+1,i,MaticeEnum)); 3118 } 3117 materials->AddObject(new Matice(i+1,i,MaticeEnum)); 3118 } 3119 3119 } 3120 3120 3121 3121 /*Add new constant material property to materials, at the end: */ 3122 3122 Matpar *newmatpar=static_cast<Matpar*>(this->materials->GetObjectByOffset(this->materials->Size()-1)->copy()); 3123 3123 newmatpar->SetMid(newnumberofelements+1); 3124 materials->AddObject(newmatpar);//put it at the end of the materials 3124 materials->AddObject(newmatpar);//put it at the end of the materials 3125 3125 } 3126 3126 /*}}}*/ 3127 3127 void FemModel::CreateNodes(int newnumberofvertices,int* my_vertices,int nodecounter,int analysis_enum,Nodes* nodes){/*{{{*/ … … 3128 3128 3129 3129 int lid=0; 3130 3130 for(int j=0;j<newnumberofvertices;j++){ 3131 if(my_vertices[j]){ 3132 3133 Node* newnode=new Node(); 3134 3131 if(my_vertices[j]){ 3132 3133 Node* newnode=new Node(); 3134 3135 3135 /*id: */ 3136 3136 newnode->id=nodecounter+j+1; 3137 3137 newnode->sid=j; 3138 3138 newnode->lid=lid++; 3139 3139 newnode->analysis_enum=analysis_enum; 3140 3140 3141 3141 /*Initialize coord_system: Identity matrix by default*/ 3142 3142 for(int k=0;k<3;k++) for(int l=0;l<3;l++) newnode->coord_system[k][l]=0.0; 3143 3143 for(int k=0;k<3;k++) newnode->coord_system[k][k]=1.0; 3144 3144 3145 3145 /*indexing:*/ 3146 3146 newnode->indexingupdate=true; 3147 3147 3148 3148 Analysis* analysis=EnumToAnalysis(analysis_enum); 3149 3149 int *doftypes=NULL; 3150 3150 int numdofs=analysis->DofsPerNode(&doftypes,Domain2DhorizontalEnum,SSAApproximationEnum); … … 3158 3158 3159 3159 /*Stressbalance Horiz*/ 3160 3160 if(analysis_enum==StressbalanceAnalysisEnum){ 3161 // itapopo this code is rarely used. 3161 // itapopo this code is rarely used. 3162 3162 /*Coordinate system provided, convert to coord_system matrix*/ 3163 3163 //XZvectorsToCoordinateSystem(&this->coord_system[0][0],&iomodel->Data(StressbalanceReferentialEnum)[j*6]); 3164 3164 //_assert_(sqrt( coord_system[0][0]*coord_system[0][0] + coord_system[1][0]*coord_system[1][0]) >1.e-4); … … 3173 3173 3174 3174 if(!femmodel_vertices) _error_("GetMesh: vertices are NULL."); 3175 3175 if(!femmodel_elements) _error_("GetMesh: elements are NULL."); 3176 3176 3177 3177 int numberofvertices = femmodel_vertices->NumberOfVertices(); 3178 3178 int numberofelements = femmodel_elements->NumberOfElements(); 3179 3179 int elementswidth = this->GetElementsWidth(); // just 2D mesh in this version (just tria elements) 3180 3180 IssmDouble* x = NULL; 3181 3181 IssmDouble* y = NULL; 3182 IssmDouble* z = NULL; 3182 IssmDouble* z = NULL; 3183 3183 int* elementslist = NULL; 3184 3184 int* elem_vertices = NULL; 3185 3185 IssmDouble *id1 = NULL; … … 3188 3188 3189 3189 /*Get vertices coordinates*/ 3190 3190 VertexCoordinatesx(&x,&y,&z,femmodel_vertices,false) ; 3191 3191 3192 3192 /*Get element vertices*/ 3193 3193 elem_vertices = xNew<int>(elementswidth); 3194 3194 Vector<IssmDouble>* vid1= new Vector<IssmDouble>(numberofelements); … … 3203 3203 vid2->SetValue(element->sid,elem_vertices[1],INS_VAL); 3204 3204 vid3->SetValue(element->sid,elem_vertices[2],INS_VAL); 3205 3205 } 3206 3206 3207 3207 /*Assemble*/ 3208 3208 vid1->Assemble(); 3209 3209 vid2->Assemble(); … … 3213 3213 id1 = vid1->ToMPISerial(); 3214 3214 id2 = vid2->ToMPISerial(); 3215 3215 id3 = vid3->ToMPISerial(); 3216 3216 3217 3217 /*Construct elements list*/ 3218 3218 elementslist=xNew<int>(numberofelements*elementswidth); 3219 3219 if(numberofelements*elementswidth<0) _error_("numberofelements negative."); … … 3222 3222 elementslist[elementswidth*i+1] = reCast<int>(id2[i])+1; //InterpMesh wants Matlab indexing 3223 3223 elementslist[elementswidth*i+2] = reCast<int>(id3[i])+1; //InterpMesh wants Matlab indexing 3224 3224 } 3225 3225 3226 3226 /*Assign pointers*/ 3227 3227 *px = x; 3228 3228 *py = y; … … 3243 3243 3244 3244 if(!femmodel_vertices) _error_("GetMesh: vertices are NULL."); 3245 3245 if(!femmodel_elements) _error_("GetMesh: elements are NULL."); 3246 3246 3247 3247 int numberofvertices = femmodel_vertices->Size(); //number of vertices of this partition 3248 3248 int numbertotalofvertices = femmodel_vertices->NumberOfVertices(); //number total of vertices (entire mesh) 3249 3249 int numberofelements = femmodel_elements->Size(); //number of elements of this partition … … 3250 3250 int elementswidth = this->GetElementsWidth(); //just 2D mesh in this version (just tria elements) 3251 3251 IssmDouble* x = NULL; 3252 3252 IssmDouble* y = NULL; 3253 IssmDouble* z = NULL; 3253 IssmDouble* z = NULL; 3254 3254 int* elementslist = NULL; 3255 3255 int* sidtoindex = NULL; 3256 3256 int* elem_vertices = NULL; 3257 3257 3258 3258 /*Get vertices coordinates of this partition*/ 3259 3259 sidtoindex = xNewZeroInit<int>(numbertotalofvertices);//entire mesh, all vertices 3260 3260 x = xNew<IssmDouble>(numberofvertices);//just this partition 3261 3261 y = xNew<IssmDouble>(numberofvertices);//just this partitio; 3262 3262 z = xNew<IssmDouble>(numberofvertices);//just this partitio; 3263 3263 3264 3264 /*Go through in this partition (vertices)*/ 3265 3265 for(int i=0;i<numberofvertices;i++){//just this partition 3266 Vertex* vertex=(Vertex*)femmodel_vertices->GetObjectByOffset(i); 3266 Vertex* vertex=(Vertex*)femmodel_vertices->GetObjectByOffset(i); 3267 3267 /*Attention: no spherical coordinates*/ 3268 3268 x[i]=vertex->GetX(); 3269 3269 y[i]=vertex->GetY(); … … 3276 3276 elem_vertices= xNew<int>(elementswidth); 3277 3277 elementslist = xNew<int>(numberofelements*elementswidth); 3278 3278 if(numberofelements*elementswidth<0) _error_("numberofelements negative."); 3279 3279 3280 3280 for(int i=0;i<numberofelements;i++){//just this partition 3281 3281 Element* element=xDynamicCast<Element*>(femmodel_elements->GetObjectByOffset(i)); 3282 3282 element->GetVerticesSidList(elem_vertices); … … 3283 3283 elementslist[elementswidth*i+0] = sidtoindex[elem_vertices[0]]+1; //InterpMesh wants Matlab indexing 3284 3284 elementslist[elementswidth*i+1] = sidtoindex[elem_vertices[1]]+1; //InterpMesh wants Matlab indexing 3285 3285 elementslist[elementswidth*i+2] = sidtoindex[elem_vertices[2]]+1; //InterpMesh wants Matlab indexing 3286 } 3287 3286 } 3287 3288 3288 /*Assign pointers*/ 3289 3289 *px = x; 3290 3290 *py = y; 3291 3291 *pz = z; 3292 3292 *pelementslist = elementslist; //Matlab indexing. InterMesh uses this type. 3293 *psidtoindex = sidtoindex; //it is ncessary to insert inputs 3293 *psidtoindex = sidtoindex; //it is ncessary to insert inputs 3294 3294 3295 3295 /*Cleanup*/ 3296 3296 xDelete<int>(elem_vertices); … … 3301 3301 /*ATTENTION: JUST SPCVX AND SPCVY*/ 3302 3302 /*OTHERS CONSTRAINTS MUST BE IMPLEMENTED*/ 3303 3303 if(analysis_enum!=StressbalanceAnalysisEnum) return; 3304 3304 3305 3305 int numberofnodes_analysistype= this->nodes->NumberOfNodes(analysis_enum); 3306 int dofpernode = 2; //vx and vy 3306 int dofpernode = 2; //vx and vy 3307 3307 int numberofcols = dofpernode*2; //to keep dofs and flags in the vspc vector 3308 3308 int numberofvertices = this->vertices->NumberOfVertices(); //global, entire old mesh 3309 3309 int numberofelements = this->elements->NumberOfElements(); //global, entire old mesh … … 3327 3327 newx=xNew<IssmDouble>(newnumberofvertices);//just the new partition 3328 3328 newy=xNew<IssmDouble>(newnumberofvertices);//just the new partition 3329 3329 for(int i=0;i<newnumberofvertices;i++){//just the new partition 3330 Vertex* vertex=(Vertex*)newfemmodel_vertices->GetObjectByOffset(i); 3330 Vertex* vertex=(Vertex*)newfemmodel_vertices->GetObjectByOffset(i); 3331 3331 /*Attention: no spherical coordinates*/ 3332 3332 newx[i]=vertex->GetX(); 3333 3333 newy[i]=vertex->GetY(); … … 3335 3335 3336 3336 /*Get spcvx and spcvy of old mesh*/ 3337 3337 for(int i=0;i<this->constraints->Size();i++){ 3338 3338 3339 3339 Constraint* constraint=(Constraint*)constraints->GetObjectByOffset(i); 3340 3340 if(!constraint->InAnalysis(analysis_enum)) _error_("AMR create constraints for "<<EnumToStringx(analysis_enum)<<" not supported yet!\n"); 3341 3341 … … 3342 3342 SpcStatic* spcstatic = xDynamicCast<SpcStatic*>(constraint); 3343 3343 int dof = spcstatic->GetDof(); 3344 3344 int node = spcstatic->GetNodeId(); 3345 IssmDouble spcvalue = spcstatic->GetValue(); 3345 IssmDouble spcvalue = spcstatic->GetValue(); 3346 3346 int nodeindex = node-1; 3347 3347 3348 3348 /*vx and vx flag insertion*/ 3349 3349 if(dof==0) {//vx 3350 3350 vspc->SetValue(nodeindex*numberofcols,spcvalue,INS_VAL); //vx 3351 3351 vspc->SetValue(nodeindex*numberofcols+dofpernode,1,INS_VAL);//vxflag 3352 } 3352 } 3353 3353 /*vy and vy flag insertion*/ 3354 3354 if(dof==1){//vy 3355 3355 vspc->SetValue(nodeindex*numberofcols+1,spcvalue,INS_VAL); //vy … … 3365 3365 InterpFromMeshToMesh2dx(&newspc,elementslist,x,y,numberofvertices,numberofelements, 3366 3366 spc,numberofvertices,numberofcols, 3367 3367 newx,newy,newnumberofvertices,NULL); 3368 3368 3369 3369 /*Now, insert the interpolated constraints in the data set (constraints)*/ 3370 3370 count=0; 3371 3371 for(int i=0;i<newnumberofvertices;i++){//just in the new partition … … 3382 3382 Vertex* vertex=(Vertex*)newfemmodel_vertices->GetObjectByOffset(i); 3383 3383 /*spcvy*/ 3384 3384 if(!xIsNan<IssmDouble>(newspc[i*numberofcols+1]) && newspc[i*numberofcols+dofpernode+1]>(1-eps) ){ 3385 newfemmodel_constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+vertex->Sid()+1,1,newspc[i*numberofcols+1],analysis_enum)); 3385 newfemmodel_constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+vertex->Sid()+1,1,newspc[i*numberofcols+1],analysis_enum)); 3386 3386 //add count'th spc, on node i+1, setting dof 1 to vx. 3387 3387 count++; 3388 3388 } … … 3439 3439 int numprocs = IssmComm::GetSize(); 3440 3440 bool *my_elements = NULL; 3441 3441 int *my_vertices = NULL; 3442 3443 _assert_(newnumberofvertices>0); 3444 _assert_(newnumberofelements>0); 3442 3443 _assert_(newnumberofvertices>0); 3444 _assert_(newnumberofelements>0); 3445 3445 epart=xNew<int>(newnumberofelements); 3446 3446 npart=xNew<int>(newnumberofvertices); 3447 3447 index=xNew<int>(elementswidth*newnumberofelements); 3448 3448 3449 3449 for (int i=0;i<newnumberofelements;i++){ 3450 3450 for (int j=0;j<elementswidth;j++){ 3451 3451 *(index+elementswidth*i+j)=(*(newelementslist+elementswidth*i+j))-1; //-1 for C indexing in Metis … … 3465 3465 for (int i=0;i<newnumberofelements;i++) epart[i]=0; 3466 3466 for (int i=0;i<newnumberofvertices;i++) npart[i]=0; 3467 3467 } 3468 else _error_("At least one processor is required"); 3468 else _error_("At least one processor is required"); 3469 3469 3470 3470 my_vertices=xNew<int>(newnumberofvertices); 3471 3471 my_elements=xNew<bool>(newnumberofelements); … … 3475 3475 /*Start figuring out, out of the partition, which elements belong to this cpu: */ 3476 3476 for(int i=0;i<newnumberofelements;i++){ 3477 3477 /*!All elements have been partitioned above, only deal with elements for this cpu: */ 3478 if(my_rank==epart[i]){ 3478 if(my_rank==epart[i]){ 3479 3479 my_elements[i]=true; 3480 /*Now that we are here, we can also start building the list of vertices belonging to this cpu partition: we use 3481 *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) 3482 into the vertices coordinates. If we start plugging 1 into my_vertices for each index[n][i] (i=0:2), then my_vertices 3480 /*Now that we are here, we can also start building the list of vertices belonging to this cpu partition: we use 3481 *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) 3482 into the vertices coordinates. If we start plugging 1 into my_vertices for each index[n][i] (i=0:2), then my_vertices 3483 3483 will hold which vertices belong to this partition*/ 3484 3484 for(int j=0;j<elementswidth;j++){ 3485 _assert_(newelementslist[elementswidth*i+j]-1<newnumberofvertices);//newelementslist is in Matlab indexing 3485 _assert_(newelementslist[elementswidth*i+j]-1<newnumberofvertices);//newelementslist is in Matlab indexing 3486 3486 my_vertices[newelementslist[elementswidth*i+j]-1]=1;//newelementslist is in Matlab indexing 3487 3487 } 3488 3488 } … … 3494 3494 3495 3495 /*Free ressources:*/ 3496 3496 xDelete<int>(epart); 3497 xDelete<int>(npart); 3497 xDelete<int>(npart); 3498 3498 xDelete<int>(index); 3499 3499 } 3500 3500 /*}}}*/ 3501 3501 void FemModel::SmoothedDeviatoricStressTensor(IssmDouble** ptauxx,IssmDouble** ptauyy,IssmDouble** ptauxy){/*{{{*/ 3502 3502 3503 3503 int elementswidth = this->GetElementsWidth();//just 2D mesh, tria elements 3504 3504 int numberofvertices = this->vertices->NumberOfVertices(); 3505 3505 IssmDouble weight = 0.; 3506 IssmDouble* tauxx = NULL; 3507 IssmDouble* tauyy = NULL; 3508 IssmDouble* tauxy = NULL; 3506 IssmDouble* tauxx = NULL; 3507 IssmDouble* tauyy = NULL; 3508 IssmDouble* tauxy = NULL; 3509 3509 IssmDouble* totalweight = NULL; 3510 3510 IssmDouble* deviatoricstressxx = xNew<IssmDouble>(elementswidth); 3511 3511 IssmDouble* deviatoricstressyy = xNew<IssmDouble>(elementswidth); … … 3515 3515 Vector<IssmDouble>* vectauyy = new Vector<IssmDouble>(numberofvertices); 3516 3516 Vector<IssmDouble>* vectauxy = new Vector<IssmDouble>(numberofvertices); 3517 3517 Vector<IssmDouble>* vectotalweight = new Vector<IssmDouble>(numberofvertices); 3518 3518 3519 3519 /*Update the Deviatoric Stress tensor over the elements*/ 3520 3520 this->DeviatoricStressx(); 3521 3521 3522 3522 /*Calculate the Smoothed Deviatoric Stress tensor*/ 3523 3523 for(int i=0;i<this->elements->Size();i++){ 3524 3524 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i)); … … 3563 3563 3564 3564 /*Divide for the total weight*/ 3565 3565 for(int i=0;i<numberofvertices;i++){ 3566 _assert_(totalweight[i]>0); 3566 _assert_(totalweight[i]>0); 3567 3567 tauxx[i] = tauxx[i]/totalweight[i]; 3568 3568 tauyy[i] = tauyy[i]/totalweight[i]; 3569 3569 tauxy[i] = tauxy[i]/totalweight[i]; … … 3588 3588 /*}}}*/ 3589 3589 void FemModel::ZZErrorEstimator(IssmDouble** pelementerror){/*{{{*/ 3590 3590 3591 /*Compute the Zienkiewicz and Zhu (ZZ) error estimator for the deviatoric stress tensor. 3591 /*Compute the Zienkiewicz and Zhu (ZZ) error estimator for the deviatoric stress tensor. 3592 3592 * Ref.: Zienkiewicz and Zhu, A Simple Error Estimator and Adaptive Procedure for Practical Engineering Analysis, Int. J. Numer. Meth. Eng, 1987*/ 3593 3593 3594 3594 IssmDouble Jdet,error,ftxx,ftyy,ftxy; … … 3608 3608 3609 3609 /*Get smoothed deviatoric stress tensor*/ 3610 3610 this->SmoothedDeviatoricStressTensor(&smoothedtauxx,&smoothedtauyy,&smoothedtauxy); 3611 3611 3612 3612 /*Integrate the error over elements*/ 3613 3613 for(int i=0;i<this->elements->Size();i++){ 3614 3614 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i)); … … 3616 3616 element->GetInputListOnVertices(tauyy,DeviatoricStressyyEnum); 3617 3617 element->GetInputListOnVertices(tauxy,DeviatoricStressxyEnum); 3618 3618 element->GetVerticesSidList(elem_vertices); 3619 3619 3620 3620 /*Integrate*/ 3621 3621 element->GetVerticesCoordinates(&xyz_list); 3622 3622 Gauss* gauss=element->NewGauss(2); … … 3631 3631 ftyy+=(tauyy[n]-smoothedtauyy[elem_vertices[n]])*basis[n]; 3632 3632 ftxy+=(tauxy[n]-smoothedtauxy[elem_vertices[n]])*basis[n]; 3633 3633 } 3634 error+=Jdet*gauss->weight*( pow(ftxx,2)+pow(ftyy,2)+pow(ftxy,2) ); //e^2 3634 error+=Jdet*gauss->weight*( pow(ftxx,2)+pow(ftyy,2)+pow(ftxy,2) ); //e^2 3635 3635 } 3636 /*Set the error in the global vector*/ 3636 /*Set the error in the global vector*/ 3637 3637 sid=element->Sid(); 3638 3638 error = sqrt(error);//sqrt(e^2) 3639 3639 velementerror->SetValue(sid,error,INS_VAL); … … 3647 3647 3648 3648 /*Serialize and set output*/ 3649 3649 (*pelementerror)=velementerror->ToMPISerial(); 3650 3650 3651 3651 /*Cleanup*/ 3652 3652 xDelete<IssmDouble>(smoothedtauxx); 3653 3653 xDelete<IssmDouble>(smoothedtauyy); … … 3691 3691 /*weight to calculate the smoothed grad H*/ 3692 3692 Tria* triaelement = xDynamicCast<Tria*>(element); 3693 3693 weight = triaelement->GetArea();//the tria area is a choice for the weight 3694 3694 3695 3695 /*dH/dx*/ 3696 3696 vecdHdx->SetValue(elem_vertices[0],weight*GradH[0],ADD_VAL); 3697 3697 vecdHdx->SetValue(elem_vertices[1],weight*GradH[0],ADD_VAL); … … 3759 3759 3760 3760 /*Get smoothed deviatoric stress tensor*/ 3761 3761 this->SmoothedGradThickness(&smoothed_dHdx,&smoothed_dHdy); 3762 3762 3763 3763 /*Integrate the error over elements*/ 3764 3764 for(int i=0;i<this->elements->Size();i++){ 3765 3765 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i)); … … 3845 3845 Vector<IssmDouble>* vyc = new Vector<IssmDouble>(numberofelements); 3846 3846 IssmDouble* x = NULL; 3847 3847 IssmDouble* y = NULL; 3848 IssmDouble* z = NULL; 3848 IssmDouble* z = NULL; 3849 3849 IssmDouble* xyz_list = NULL; 3850 3850 IssmDouble x1,y1,x2,y2,x3,y3; 3851 3851 … … 3854 3854 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i)); 3855 3855 //element->GetVerticesSidList(elem_vertices); 3856 3856 int sid = element->Sid(); 3857 element->GetVerticesCoordinates(&xyz_list); 3857 element->GetVerticesCoordinates(&xyz_list); 3858 3858 x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1]; 3859 3859 x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1]; 3860 3860 x3 = xyz_list[3*2+0];y3 = xyz_list[3*2+1]; … … 3887 3887 if(levelset_type!=MaskGroundediceLevelsetEnum && levelset_type!=MaskIceLevelsetEnum){ 3888 3888 _error_("level set type not implemented yet!"); 3889 3889 } 3890 3890 3891 3891 /*Outputs*/ 3892 3892 IssmDouble* zerolevelset_points = NULL; 3893 3893 int npoints = 0; 3894 3894 3895 3895 /*Intermediaries*/ 3896 3896 int elementswidth = this->GetElementsWidth(); 3897 3897 int numberofelements = this->elements->NumberOfElements(); … … 3904 3904 IssmDouble* y_zerolevelset = NULL; 3905 3905 int count,sid; 3906 3906 IssmDouble xc,yc,x1,y1,x2,y2,x3,y3; 3907 3907 3908 3908 /*Use the element center coordinate if level set is zero (grounding line or ice front), otherwise set NAN*/ 3909 3909 for(int i=0;i<this->elements->Size();i++){ 3910 3910 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i)); … … 3911 3911 element->GetInputListOnVertices(levelset,levelset_type); 3912 3912 element->GetVerticesSidList(elem_vertices); 3913 3913 sid= element->Sid(); 3914 element->GetVerticesCoordinates(&xyz_list); 3914 element->GetVerticesCoordinates(&xyz_list); 3915 3915 x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1]; 3916 3916 x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1]; 3917 3917 x3 = xyz_list[3*2+0];y3 = xyz_list[3*2+1]; 3918 3918 xc = NAN; 3919 yc = NAN; 3919 yc = NAN; 3920 3920 Tria* tria = xDynamicCast<Tria*>(element); 3921 3921 if(tria->IsIceInElement()){/*verify if there is ice in the element*/ 3922 if(levelset[0]*levelset[1]<0. || levelset[0]*levelset[2]<0. || 3922 if(levelset[0]*levelset[1]<0. || levelset[0]*levelset[2]<0. || 3923 3923 abs(levelset[0]*levelset[1])<DBL_EPSILON || abs(levelset[0]*levelset[2])<DBL_EPSILON) { 3924 3924 xc=(x1+x2+x3)/3.; 3925 3925 yc=(y1+y2+y3)/3.; … … 3949 3949 count++; 3950 3950 } 3951 3951 } 3952 3952 3953 3953 /*Assign outputs*/ 3954 3954 numberofpoints = npoints; 3955 3955 (*pzerolevelset_points) = zerolevelset_points; … … 3989 3989 /*save the d_responses pointer: */ 3990 3990 responses_pointer=d_responses; 3991 3991 3992 //watch out, we have more d_numresponses than numresponsedescriptors, because the responses have been expanded if they were scaled. 3992 //watch out, we have more d_numresponses than numresponsedescriptors, because the responses have been expanded if they were scaled. 3993 3993 //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: */ 3994 3994 3995 3995 for(i=0;i<numresponsedescriptors;i++){ … … 4082 4082 void FemModel::EsaGeodetic2D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pX, Vector<IssmDouble>* pY, IssmDouble* xx, IssmDouble* yy){/*{{{*/ 4083 4083 4084 4084 int ns,nsmax; 4085 4085 4086 4086 /*Go through elements, and add contribution from each element to the deflection vector wg:*/ 4087 4087 ns = elements->Size(); 4088 4088 4089 4089 /*Figure out max of ns: */ 4090 4090 ISSM_MPI_Reduce(&ns,&nsmax,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm()); 4091 4091 ISSM_MPI_Bcast(&nsmax,1,ISSM_MPI_INT,0,IssmComm::GetComm()); … … 4104 4104 pY->Assemble(); 4105 4105 } 4106 4106 } 4107 4107 4108 4108 /*One last time: */ 4109 4109 pUp->Assemble(); 4110 4110 pNorth->Assemble(); … … 4123 4123 IssmDouble eartharea_cpu=0; 4124 4124 4125 4125 int ns,nsmax; 4126 4126 4127 4127 /*Go through elements, and add contribution from each element to the deflection vector wg:*/ 4128 4128 ns = elements->Size(); 4129 4130 /*First, figure out the surface area of Earth: */ 4129 4130 /*First, figure out the surface area of Earth: */ 4131 4131 for(int i=0;i<ns;i++){ 4132 4132 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i)); 4133 4133 eartharea_cpu += element->GetAreaSpherical(); … … 4151 4151 pEast->Assemble(); 4152 4152 } 4153 4153 } 4154 4154 4155 4155 /*One last time: */ 4156 4156 pUp->Assemble(); 4157 4157 pNorth->Assemble(); … … 4213 4213 if(i%loop==0)pRSLgi->Assemble(); 4214 4214 } 4215 4215 if(VerboseConvergence())_printf0_("\n"); 4216 4216 4217 4217 /*One last time: */ 4218 4218 pRSLgi->Assemble(); 4219 4219 … … 4231 4231 4232 4232 /*serialized vectors:*/ 4233 4233 IssmDouble* RSLg_old=NULL; 4234 4234 4235 4235 IssmDouble eartharea=0; 4236 4236 IssmDouble eartharea_cpu=0; 4237 4237 4238 4238 int ns,nsmax; 4239 4239 4240 4240 /*Serialize vectors from previous iteration:*/ 4241 4241 RSLg_old=pRSLg_old->ToMPISerial(); 4242 4242 … … 4248 4248 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i)); 4249 4249 eartharea_cpu += element->GetAreaSpherical(); 4250 4250 } 4251 4251 4252 4252 ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); 4253 4253 ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 4254 4254 … … 4266 4266 if(i%loop==0)pRSLgo->Assemble(); 4267 4267 } 4268 4268 if(verboseconvolution)if(VerboseConvergence())_printf0_("\n"); 4269 4269 4270 4270 /*Free ressources:*/ 4271 4271 xDelete<IssmDouble>(RSLg_old); 4272 4272 } … … 4354 4354 4355 4355 /*Free ressources:*/ 4356 4356 xDelete<IssmDouble>(RSLg_old); 4357 4358 4357 4358 4359 4359 } 4360 4360 /*}}}*/ 4361 4361 void FemModel::SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pRSLg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/ … … 4362 4362 4363 4363 /*serialized vectors:*/ 4364 4364 IssmDouble* RSLg=NULL; 4365 4365 4366 4366 IssmDouble eartharea=0; 4367 4367 IssmDouble eartharea_cpu=0; 4368 4368 4369 4369 int ns,nsmax; 4370 4370 4371 4371 /*Serialize vectors from previous iteration:*/ 4372 4372 RSLg=pRSLg->ToMPISerial(); 4373 4373 4374 4374 /*Go through elements, and add contribution from each element to the deflection vector wg:*/ 4375 4375 ns = elements->Size(); 4376 4376 4377 4377 /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */ 4378 4378 for(int i=0;i<ns;i++){ 4379 4379 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i)); … … 4401 4401 } 4402 4402 } 4403 4403 } 4404 4404 4405 4405 /*One last time: */ 4406 4406 pUp->Assemble(); 4407 4407 if(horiz){ … … 4435 4435 } 4436 4436 ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); 4437 4437 ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 4438 4438 4439 4439 ISSM_MPI_Reduce (&oceanvalue_cpu,&oceanvalue,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); 4440 4440 ISSM_MPI_Bcast(&oceanvalue,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 4441 4441 4442 4442 /*Free ressources:*/ 4443 4443 xDelete<IssmDouble>(RSLg_serial); 4444 4444 4445 4445 return oceanvalue/oceanarea; 4446 4446 } 4447 4447 /*}}}*/ … … 4457 4457 IssmDouble* old_active = NULL; 4458 4458 int* eplzigzag_counter = NULL; 4459 4459 int eplflip_lock; 4460 4460 4461 4461 HydrologyDCEfficientAnalysis* effanalysis = new HydrologyDCEfficientAnalysis(); 4462 4462 HydrologyDCInefficientAnalysis* inefanalysis = new HydrologyDCInefficientAnalysis(); 4463 4463 … … 4464 4464 /*Step 1: update mask, the mask might be extended by residual and/or using downstream sediment head*/ 4465 4465 mask=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum)); 4466 4466 recurence=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum)); 4467 this->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum); 4468 this->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum); 4467 this->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum); 4468 this->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum); 4469 4469 GetVectorFromInputsx(&old_active,this,HydrologydcMaskEplactiveNodeEnum,NodeSIdEnum); 4470 4470 4471 4471 for (int i=0;i<elements->Size();i++){ 4472 4472 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i)); 4473 4473 effanalysis->HydrologyEPLGetMask(mask,recurence,eplzigzag_counter,element); … … 4485 4485 this->parameters->SetParam(eplzigzag_counter,this->nodes->Size(),EplZigZagCounterEnum); 4486 4486 /*Assemble and serialize*/ 4487 4487 mask->Assemble(); 4488 serial_mask=mask->ToMPISerial(); 4489 4488 serial_mask=mask->ToMPISerial(); 4489 4490 4490 xDelete<int>(eplzigzag_counter); 4491 4491 xDelete<IssmDouble>(serial_rec); 4492 4492 xDelete<IssmDouble>(old_active); … … 4528 4528 delete inefanalysis; 4529 4529 int sum_counter; 4530 4530 ISSM_MPI_Reduce(&counter,&sum_counter,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() ); 4531 ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 4531 ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 4532 4532 counter=sum_counter; 4533 4533 *pEplcount = counter; 4534 4534 if(VerboseSolution()) _printf0_(" Number of active nodes in EPL layer: "<< counter <<"\n"); … … 4649 4649 xDelete<IssmDouble>(serial_active); 4650 4650 int sum_counter; 4651 4651 ISSM_MPI_Reduce(&counter,&sum_counter,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() ); 4652 ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 4652 ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 4653 4653 counter=sum_counter; 4654 4654 *pL2count = counter; 4655 4655 if(VerboseSolution()) _printf0_(" Number of active nodes L2 Projection: "<< counter <<"\n"); … … 4734 4734 } 4735 4735 } 4736 4736 /*}}}*/ 4737 #ifdef _HAVE_JAVASCRIPT_ 4737 #ifdef _HAVE_JAVASCRIPT_ 4738 4738 FemModel::FemModel(IssmDouble* buffer, int buffersize, char* toolkits, char* solution, char* modelname,ISSM_MPI_Comm incomm, bool trace){ /*{{{*/ 4739 4739 /*configuration: */ 4740 4740 int solution_type; … … 4749 4749 4750 4750 /*From command line arguments, retrieve different filenames needed to create the FemModel: */ 4751 4751 solution_type=StringToEnumx(solution); 4752 4752 4753 4753 /*Create femmodel from input files: */ 4754 4754 profiler->Start(MPROCESSOR); 4755 4755 this->InitFromBuffers((char*)buffer,buffersize,toolkits, solution_type,trace,NULL); 4756 4756 profiler->Stop(MPROCESSOR); 4757 4757 4758 4758 /*Save communicator in the parameters dataset: */ 4759 4759 this->parameters->AddObject(new GenericParam<ISSM_MPI_Comm>(incomm,FemModelCommEnum)); 4760 4760 … … 4769 4769 char** poutputbuffer; 4770 4770 size_t* poutputbuffersize; 4771 4771 4772 4772 4773 4773 /*Before we delete the profiler, report statistics for this run: */ 4774 4774 profiler->Stop(TOTAL); //final tagging 4775 4775 _printf0_("\n"); … … 4782 4782 <<profiler->TotalTimeModSec(TOTAL)<<" sec" 4783 4783 ); 4784 4784 _printf0_("\n"); 4785 4785 4786 4786 /*Before we close the output file, recover the buffer and size:*/ 4787 4787 outputbufferparam = xDynamicCast<GenericParam<char**>*>(this->parameters->FindParamObject(OutputBufferPointerEnum)); 4788 4788 poutputbuffer=outputbufferparam->GetParameterValue(); … … 4822 4822 fclose(toolkitsoptionsfid); 4823 4823 4824 4824 /*Open output file once for all and add output file descriptor to parameters*/ 4825 output_fid=open_memstream(&outputbuffer,&outputsize); 4826 if(output_fid==NULL ||output_fid==0)_error_("could not initialize output stream");4825 output_fid=open_memstream(&outputbuffer,&outputsize); 4826 if(output_fid==NULL)_error_("could not initialize output stream"); 4827 4827 this->parameters->SetParam(output_fid,OutputFilePointerEnum); 4828 4828 this->parameters->AddObject(new GenericParam<char**>(&outputbuffer,OutputBufferPointerEnum)); 4829 4829 this->parameters->AddObject(new GenericParam<size_t*>(&outputsize,OutputBufferSizePointerEnum)); … … 4864 4864 this->amrbamg->thicknesserror_threshold>0||this->amrbamg->deviatoricerror_threshold>0){ 4865 4865 /*Initialize hmaxvertices with NAN*/ 4866 4866 hmaxvertices_serial=xNew<IssmDouble>(numberofvertices); 4867 for(int i=0;i<numberofvertices;i++) hmaxvertices_serial[i]=NAN; 4867 for(int i=0;i<numberofvertices;i++) hmaxvertices_serial[i]=NAN; 4868 4868 /*Fill hmaxvertices*/ 4869 4869 if(this->amrbamg->groundingline_distance>0) this->GethmaxVerticesFromZeroLevelSetDistance(hmaxvertices_serial,MaskGroundediceLevelsetEnum); 4870 4870 if(this->amrbamg->icefront_distance>0) this->GethmaxVerticesFromZeroLevelSetDistance(hmaxvertices_serial,MaskIceLevelsetEnum); … … 4871 4871 if(this->amrbamg->thicknesserror_threshold>0) this->GethmaxVerticesFromEstimators(hmaxvertices_serial,ThicknessErrorEstimatorEnum); 4872 4872 if(this->amrbamg->deviatoricerror_threshold>0) this->GethmaxVerticesFromEstimators(hmaxvertices_serial,DeviatoricStressErrorEstimatorEnum); 4873 4873 } 4874 4874 4875 4875 if(my_rank==0){ 4876 4876 this->amrbamg->ExecuteRefinementBamg(vector_serial,hmaxvertices_serial,&newnumberofvertices,&newnumberofelements,&newx,&newy,&newz,&newelementslist); 4877 4877 if(newnumberofvertices<=0 || newnumberofelements<=0) _error_("Error in the refinement process."); … … 4891 4891 newz=xNew<IssmDouble>(newnumberofvertices); 4892 4892 newelementslist=xNew<int>(newnumberofelements*this->GetElementsWidth()); 4893 4893 } 4894 ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 4895 ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 4896 ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 4897 ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm()); 4894 ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 4895 ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 4896 ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 4897 ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm()); 4898 4898 4899 4899 /*Assign output pointers*/ 4900 4900 *pnewnumberofvertices = newnumberofvertices; … … 4927 4927 4928 4928 /*Create bamg data structures for bamg*/ 4929 4929 this->amrbamg = new AmrBamg(); 4930 4930 4931 4931 /*Get amr parameters*/ 4932 4932 this->parameters->FindParam(&hmin,AmrHminEnum); 4933 4933 this->parameters->FindParam(&hmax,AmrHmaxEnum); … … 4951 4951 this->amrbamg->SetBamgOpts(hmin,hmax,err,gradation); 4952 4952 4953 4953 /*Re-create original mesh and put it in bamg structure (only cpu 0)*/ 4954 if(my_rank==0){ 4954 if(my_rank==0){ 4955 4955 this->amrbamg->Initialize(elements,x,y,numberofvertices,numberofelements); 4956 4956 } 4957 4957 … … 4965 4965 void FemModel::GethmaxVerticesFromZeroLevelSetDistance(IssmDouble* hmaxvertices,int levelset_type){/*{{{*/ 4966 4966 4967 4967 if(!hmaxvertices) _error_("hmaxvertices is NULL!\n"); 4968 4968 4969 4969 /*Intermediaries*/ 4970 4970 int numberofvertices = this->vertices->NumberOfVertices(); 4971 4971 IssmDouble* verticedistance = NULL; … … 4972 4972 IssmDouble threshold,resolution; 4973 4973 4974 4974 switch(levelset_type){ 4975 case MaskGroundediceLevelsetEnum: 4975 case MaskGroundediceLevelsetEnum: 4976 4976 threshold = this->amrbamg->groundingline_distance; 4977 4977 resolution = this->amrbamg->groundingline_resolution; 4978 4978 break; … … 4986 4986 /*Get vertice distance to zero level set points*/ 4987 4987 this->GetVerticeDistanceToZeroLevelSet(&verticedistance,levelset_type); 4988 4988 if(!verticedistance) _error_("verticedistance is NULL!\n"); 4989 4989 4990 4990 /*Fill hmaxVertices*/ 4991 4991 for(int i=0;i<numberofvertices;i++){ 4992 4992 if(verticedistance[i]<threshold){ … … 5000 5000 } 5001 5001 /*}}}*/ 5002 5002 void FemModel::GethmaxVerticesFromEstimators(IssmDouble* hmaxvertices,int errorestimator_type){/*{{{*/ 5003 5003 5004 5004 if(!hmaxvertices) _error_("hmaxvertices is NULL!\n"); 5005 5005 5006 5006 /*Intermediaries*/ … … 5008 5008 int numberofelements = this->elements->NumberOfElements(); 5009 5009 int numberofvertices = this->vertices->NumberOfVertices(); 5010 5010 IssmDouble* maxlength = xNew<IssmDouble>(numberofelements); 5011 IssmDouble* error_vertices = xNewZeroInit<IssmDouble>(numberofvertices); 5011 IssmDouble* error_vertices = xNewZeroInit<IssmDouble>(numberofvertices); 5012 5012 IssmDouble* error_elements = NULL; 5013 5013 IssmDouble* x = NULL; 5014 5014 IssmDouble* y = NULL; … … 5021 5021 5022 5022 /*Fill variables*/ 5023 5023 switch(errorestimator_type){ 5024 case ThicknessErrorEstimatorEnum: 5024 case ThicknessErrorEstimatorEnum: 5025 5025 threshold = this->amrbamg->thicknesserror_threshold; 5026 5026 groupthreshold = this->amrbamg->thicknesserror_groupthreshold; 5027 5027 resolution = this->amrbamg->thicknesserror_resolution; … … 5046 5046 switch(errorestimator_type){ 5047 5047 case ThicknessErrorEstimatorEnum: this->amrbamg->thicknesserror_maximum = maxerror;break; 5048 5048 case DeviatoricStressErrorEstimatorEnum: this->amrbamg->deviatoricerror_maximum = maxerror;break; 5049 } 5049 } 5050 5050 } 5051 5051 5052 5052 /*Get mesh*/ 5053 5053 this->GetMesh(this->vertices,this->elements,&x,&y,&z,&index); 5054 5054 5055 5055 /*Fill error_vertices (this is the sum of all elements connected to the vertex)*/ 5056 5056 for(int i=0;i<numberofelements;i++){ 5057 5057 v1=index[i*elementswidth+0]-1;//Matlab to C indexing … … 5065 5065 error_vertices[v1]+=error_elements[i]; 5066 5066 error_vertices[v2]+=error_elements[i]; 5067 5067 error_vertices[v3]+=error_elements[i]; 5068 } 5068 } 5069 5069 5070 5070 /*Fill hmaxvertices with the criteria*/ 5071 5071 for(int i=0;i<numberofelements;i++){ … … 5079 5079 if(error_vertices[vid]>groupthreshold*maxerror) refine=true; 5080 5080 } 5081 5081 } 5082 /*Now, fill the hmaxvertices if requested*/ 5082 /*Now, fill the hmaxvertices if requested*/ 5083 5083 if(refine){ 5084 5084 for(int j=0;j<elementswidth;j++){ 5085 5085 vid=index[i*elementswidth+j]-1;//Matlab to C indexing … … 5109 5109 5110 5110 /*Output*/ 5111 5111 IssmDouble* verticedistance; 5112 5112 5113 5113 /*Intermediaries*/ 5114 5114 int numberofvertices = this->vertices->NumberOfVertices(); 5115 5115 IssmDouble* levelset_points= NULL; … … 5121 5121 5122 5122 /*Get vertices coordinates*/ 5123 5123 VertexCoordinatesx(&x,&y,&z,this->vertices,false) ; 5124 5125 /*Get points which level set is zero (center of elements with zero level set)*/ 5124 5125 /*Get points which level set is zero (center of elements with zero level set)*/ 5126 5126 this->GetZeroLevelSetPoints(&levelset_points,numberofpoints,levelset_type); 5127 5127 5128 5128 /*Find the minimal vertice distance to the zero levelset (grounding line or ice front)*/ … … 5131 5131 verticedistance[i]=INFINITY; 5132 5132 for(int j=0;j<numberofpoints;j++){ 5133 5133 distance=sqrt((x[i]-levelset_points[2*j])*(x[i]-levelset_points[2*j])+(y[i]-levelset_points[2*j+1])*(y[i]-levelset_points[2*j+1])); 5134 verticedistance[i]=min(distance,verticedistance[i]); 5134 verticedistance[i]=min(distance,verticedistance[i]); 5135 5135 } 5136 } 5136 } 5137 5137 5138 5138 /*Assign the pointer*/ 5139 5139 (*pverticedistance)=verticedistance; … … 5167 5167 /*Get fields, if requested*/ 5168 5168 if(this->amr->groundingline_distance>0) this->GetElementDistanceToZeroLevelSet(&gl_distance,MaskGroundediceLevelsetEnum); 5169 5169 if(this->amr->icefront_distance>0) this->GetElementDistanceToZeroLevelSet(&if_distance,MaskIceLevelsetEnum); 5170 if(this->amr->thicknesserror_threshold>0) this->ThicknessZZErrorEstimator(&thicknesserror); 5171 if(this->amr->deviatoricerror_threshold>0) this->ZZErrorEstimator(&deviatoricerror); 5172 5170 if(this->amr->thicknesserror_threshold>0) this->ThicknessZZErrorEstimator(&thicknesserror); 5171 if(this->amr->deviatoricerror_threshold>0) this->ZZErrorEstimator(&deviatoricerror); 5172 5173 5173 if(my_rank==0){ 5174 5174 this->amr->ExecuteRefinement(gl_distance,if_distance,deviatoricerror,thicknesserror, 5175 &newnumberofvertices,&newnumberofelements,&newx,&newy,&newelementslist); 5175 &newnumberofvertices,&newnumberofelements,&newx,&newy,&newelementslist); 5176 5176 newz=xNewZeroInit<IssmDouble>(newnumberofvertices); 5177 5177 if(newnumberofvertices<=0 || newnumberofelements<=0) _error_("Error in the ReMeshNeopz."); 5178 5178 } … … 5186 5186 newz=xNew<IssmDouble>(newnumberofvertices); 5187 5187 newelementslist=xNew<int>(newnumberofelements*this->GetElementsWidth()); 5188 5188 } 5189 ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 5190 ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 5191 ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 5192 ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm()); 5189 ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 5190 ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 5191 ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 5192 ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm()); 5193 5193 5194 /*Assign the pointers*/ 5194 /*Assign the pointers*/ 5195 5195 (*pnewelementslist) = newelementslist; //Matlab indexing 5196 5196 (*pnewx) = newx; 5197 5197 (*pnewy) = newy; … … 5207 5207 } 5208 5208 /*}}}*/ 5209 5209 void FemModel::InitializeAdaptiveRefinementNeopz(void){/*{{{*/ 5210 5210 5211 5211 /*Define variables*/ 5212 5212 int my_rank = IssmComm::GetRank(); 5213 5213 int numberofvertices = this->vertices->NumberOfVertices(); … … 5216 5216 IssmDouble* y = NULL; 5217 5217 IssmDouble* z = NULL; 5218 5218 int* elements = NULL; 5219 5219 5220 5220 /*Initialize field as NULL for now*/ 5221 5221 this->amr = NULL; 5222 5222 … … 5224 5224 /*elements comes in Matlab indexing*/ 5225 5225 this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elements); 5226 5226 5227 /*Create initial mesh (coarse mesh) in neopz data structure*/ 5227 /*Create initial mesh (coarse mesh) in neopz data structure*/ 5228 5228 /*Just CPU #0 should keep AMR object*/ 5229 5229 /*Initialize refinement pattern*/ 5230 5230 this->SetRefPatterns(); 5231 5231 this->amr = new AdaptiveMeshRefinement(); 5232 this->amr->refinement_type=1;//1 is refpattern; 0 is uniform (faster) 5232 this->amr->refinement_type=1;//1 is refpattern; 0 is uniform (faster) 5233 5233 /*Get amr parameters*/ 5234 5234 this->parameters->FindParam(&this->amr->level_max,AmrLevelMaxEnum); 5235 5235 this->parameters->FindParam(&this->amr->gradation,AmrGradationEnum); … … 5242 5242 this->parameters->FindParam(&this->amr->deviatoricerror_threshold,AmrDeviatoricErrorThresholdEnum); 5243 5243 this->parameters->FindParam(&this->amr->deviatoricerror_groupthreshold,AmrDeviatoricErrorGroupThresholdEnum); 5244 5244 this->parameters->FindParam(&this->amr->deviatoricerror_maximum,AmrDeviatoricErrorMaximumEnum); 5245 if(my_rank==0){ 5245 if(my_rank==0){ 5246 5246 this->amr->CreateInitialMesh(numberofvertices,numberofelements,x,y,elements); 5247 5247 } 5248 5248 … … 5263 5263 5264 5264 /*Output*/ 5265 5265 IssmDouble* elementdistance; 5266 5266 5267 5267 /*Intermediaries*/ 5268 5268 int numberofelements = this->elements->NumberOfElements(); 5269 5269 Vector<IssmDouble>* velementdistance = new Vector<IssmDouble>(numberofelements); … … 5272 5272 IssmDouble mindistance,distance; 5273 5273 IssmDouble xc,yc,x1,y1,x2,y2,x3,y3; 5274 5274 int numberofpoints; 5275 5276 /*Get points which level set is zero (center of elements with zero level set, levelset_points is serial)*/ 5275 5276 /*Get points which level set is zero (center of elements with zero level set, levelset_points is serial)*/ 5277 5277 this->GetZeroLevelSetPoints(&levelset_points,numberofpoints,levelset_type); 5278 5278 5279 5279 for(int i=0;i<this->elements->Size();i++){//parallel 5280 5280 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i)); 5281 5281 int sid = element->Sid(); 5282 element->GetVerticesCoordinates(&xyz_list); 5282 element->GetVerticesCoordinates(&xyz_list); 5283 5283 x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1]; 5284 5284 x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1]; 5285 5285 x3 = xyz_list[3*2+0];y3 = xyz_list[3*2+1]; … … 5289 5289 /*Loop over each point (where level set is zero)*/ 5290 5290 for(int j=0;j<numberofpoints;j++){ 5291 5291 distance =sqrt((xc-levelset_points[2*j])*(xc-levelset_points[2*j])+(yc-levelset_points[2*j+1])*(yc-levelset_points[2*j+1])); 5292 mindistance=min(distance,mindistance); 5292 mindistance=min(distance,mindistance); 5293 5293 } 5294 5294 velementdistance->SetValue(sid,mindistance,INS_VAL); 5295 5295 xDelete<IssmDouble>(xyz_list); 5296 } 5296 } 5297 5297 5298 5298 /*Assemble*/ 5299 5299 velementdistance->Assemble();
Note:
See TracBrowser
for help on using the repository browser.
