Index: /issm/trunk/src/c/modules/FieldDepthAveragex/FieldDepthAveragex.cpp =================================================================== --- /issm/trunk/src/c/modules/FieldDepthAveragex/FieldDepthAveragex.cpp (revision 3904) +++ /issm/trunk/src/c/modules/FieldDepthAveragex/FieldDepthAveragex.cpp (revision 3904) @@ -0,0 +1,48 @@ +/*!\file FieldDepthAveragex + * \brief: average field throfieldh thickness + */ + +#include "./FieldDepthAveragex.h" + +#include "../shared/shared.h" +#include "../include/include.h" +#include "../toolkits/toolkits.h" +#include "../EnumDefinitions/EnumDefinitions.h" + +void FieldDepthAveragex( Vec field, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,char* fieldname){ + + double* field_serial=NULL; + + /*First, get elements and nodes configured: */ + elements->Configure(elements,loads, nodes,vertices, materials,parameters); + nodes->Configure(elements,loads, nodes,vertices, materials,parameters); + parameters->Configure(elements,loads, nodes,vertices, materials,parameters); + + /*Serialize field: */ + VecToMPISerial(&field_serial,field); + + /*Depth average field, onto base of mesh: */ + nodes->FieldDepthAverageAtBase(field,field_serial,fieldname); + + /*Assemble vector: */ + VecAssemblyBegin(field); + VecAssemblyEnd(field); + + /*ok, field is now correct at the base. We need to extrude the value of the base + * through thickness :*/ + + /*Serialize field: */ + xfree((void**)&field_serial); + VecToMPISerial(&field_serial,field); + + /*Extrude field to all nodes: */ + nodes->FieldExtrude(field,field_serial,fieldname,0); + + /*Assemble vector: */ + VecAssemblyBegin(field); + VecAssemblyEnd(field); + + /*Free ressources:*/ + xfree((void**)&field_serial); + +} Index: /issm/trunk/src/c/modules/FieldDepthAveragex/FieldDepthAveragex.h =================================================================== --- /issm/trunk/src/c/modules/FieldDepthAveragex/FieldDepthAveragex.h (revision 3904) +++ /issm/trunk/src/c/modules/FieldDepthAveragex/FieldDepthAveragex.h (revision 3904) @@ -0,0 +1,14 @@ +/*!\file: FieldDepthAveragex.h + * \brief header file for averaging a field throfieldh thickness + */ + +#ifndef _FIELDDEPTHAVERAGEX_H +#define _FIELDDEPTHAVERAGEX_H + +#include "../DataSet/DataSet.h" + +/* local prototypes: */ +void FieldDepthAveragex( Vec field, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,char* fieldname); + +#endif /* _FIELDDEPTHAVERAGEX_H */ + Index: /issm/trunk/src/c/modules/GetSolutionFromInputsx/GetSolutionFromInputsx.cpp =================================================================== --- /issm/trunk/src/c/modules/GetSolutionFromInputsx/GetSolutionFromInputsx.cpp (revision 3904) +++ /issm/trunk/src/c/modules/GetSolutionFromInputsx/GetSolutionFromInputsx.cpp (revision 3904) @@ -0,0 +1,44 @@ +/*!\file GetSolutionFromInputsx + * \brief: update datasets using parameter inputs + */ + +#include "./GetSolutionFromInputsx.h" +#include "../shared/shared.h" +#include "../include/include.h" +#include "../toolkits/toolkits.h" +#include "../EnumDefinitions/EnumDefinitions.h" + +void GetSolutionFromInputsx( Vec* psolution, DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials, Parameters* parameters, int analysis_type, int sub_analysis_type){ + + /*intermediary: */ + int gsize; + + /*output: */ + Vec solution=NULL; + + /*First, get elements and loads configured: */ + elements->Configure(elements,loads, nodes,vertices, materials,parameters); + loads->Configure(elements, loads, nodes,vertices, materials,parameters); + nodes->Configure(elements, loads, nodes,vertices, materials,parameters); + parameters->Configure(elements,loads, nodes,vertices, materials,parameters); + + /*Get size of vector: */ + gsize=nodes->NumberOfDofs(); + + /*Initialize solution: */ + solution=NewVec(gsize); + + /*Go through elements and plug solution: */ + elements->GetSolutionFromInputs(solution,analysis_type,sub_analysis_type); + + /*Elements drive the update: */ + elements->GetSolutionFromInputs(solution,analysis_type,sub_analysis_type); + + /*Assemble vector: */ + VecAssemblyBegin(solution); + VecAssemblyEnd(solution); + + /*Assign output pointers:*/ + *psolution=solution; + +} Index: /issm/trunk/src/c/modules/GetSolutionFromInputsx/GetSolutionFromInputsx.h =================================================================== --- /issm/trunk/src/c/modules/GetSolutionFromInputsx/GetSolutionFromInputsx.h (revision 3904) +++ /issm/trunk/src/c/modules/GetSolutionFromInputsx/GetSolutionFromInputsx.h (revision 3904) @@ -0,0 +1,15 @@ +/*!\file: GetSolutionFromInputsx.h + * \brief header file for updating datasets from inputs + */ + +#ifndef _GETSOLUTIONFROMINPUTSXX_H +#define _GETSOLUTIONFROMINPUTSXX_H + +#include "../objects/objects.h" +#include "../DataSet/DataSet.h" + +/* local prototypes: */ +void GetSolutionFromInputsx( Vec* psolution, DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials, Parameters* parameters, int analysis_type, int sub_analysis_type); + +#endif /* _GETSOLUTIONFROMINPUTSXX_H */ + Index: /issm/trunk/src/c/modules/InterpFromMeshToGridx/InterpFromMeshToGridx.cpp =================================================================== --- /issm/trunk/src/c/modules/InterpFromMeshToGridx/InterpFromMeshToGridx.cpp (revision 3904) +++ /issm/trunk/src/c/modules/InterpFromMeshToGridx/InterpFromMeshToGridx.cpp (revision 3904) @@ -0,0 +1,185 @@ +/*!\file: InterpFromMeshToGridx.cpp + * \brief "c" core code for interpolating values from a structured grid. + */ + +#include "./InterpFromMeshToGridx.h" +#include "../shared/shared.h" +#include "../include/include.h" + +void InterpFromMeshToGridx(double** px_m,double** py_m,double** pgriddata,double* index_mesh, double* x_mesh, double* y_mesh, int nods,int nels, double* data_mesh, int data_length, double cornereast,double cornernorth,double xposting,double yposting,int nlines,int ncols,double default_value) { + + /*Output*/ + double* griddata=NULL; + double* x_m=NULL; + double* y_m=NULL; + + /*Intermediary*/ + int i,j,n; + int i1,i2,j1,j2; + int interpolation_type; + bool debug; + int xflip,yflip; + double area; + double area_1,area_2,area_3; + double x_tria_min,y_tria_min; + double x_tria_max,y_tria_max; + double x_grid_min,y_grid_min; + double x_grid_max,y_grid_max; + double data_value; + double cornersouth,cornerwest; + double* x_grid=NULL; + double* y_grid=NULL; + + /*some checks*/ + if (nels<1 || nods<3 || nlines<1 || ncols<1 || xposting==0 || yposting==0){ + ISSMERROR("nothing to be done according to the mesh given in input"); + } + + /*figure out what kind of interpolation is needed*/ + if (data_length==nods){ + interpolation_type=1; + } + else if (data_length==nels){ + interpolation_type=2; + } + else{ + ISSMERROR("length of vector data not supported yet. It should be of length (number of nodes) or (number of elements)!"); + } + + /*First, allocate pointers: */ + griddata=(double*)xcalloc(nlines*ncols,sizeof(double)); + x_m=(double*)xcalloc(ncols+1,sizeof(double)); + y_m=(double*)xcalloc(nlines+1,sizeof(double)); + x_grid=(double*)xcalloc(ncols,sizeof(double)); + y_grid=(double*)xcalloc(nlines,sizeof(double)); + + /*Set debug to 1 if there are lots of elements*/ + debug=(bool)((double)ncols*nlines*nels >= pow((double)10,(double)10)); + + /*figure out if x or y are flipped*/ + if (xposting<0) xflip=1; + else xflip=0; + if (yposting<0) yflip=1; + else yflip=0; + + /*Compute coordinates lists*/ + // EAST = ACTUAL WEST !!!!!!!!!!!!!!! + cornersouth=cornernorth-nlines*yposting; + cornerwest =cornereast +ncols *xposting; + for(i=0;ix_tria_max) x_tria_max=x_mesh[(int)index_mesh[3*n+i]-1]; + if(y_mesh[(int)index_mesh[3*n+i]-1]y_tria_max) y_tria_max=y_mesh[(int)index_mesh[3*n+i]-1]; + } + + /*if the current triangle is not in the grid, continue*/ + if ( (x_tria_min>x_grid_max) || (x_tria_maxy_grid_max) || (y_tria_maxy_tria_max) || (y_grid[i]x_tria_max) || (x_grid[j]=0 && area_2>=0 && area_3>=0){ + + /*Yes ! compute the value on the point*/ + if (interpolation_type==1){ + /*nodal interpolation*/ + data_value=area_1*data_mesh[(int)index_mesh[3*n+0]-1]+area_2*data_mesh[(int)index_mesh[3*n+1]-1]+area_3*data_mesh[(int)index_mesh[3*n+2]-1]; + } + else{ + /*element interpolation*/ + data_value=data_mesh[n]; + } + if (isnan(data_value)) data_value=default_value; + + /*insert value and go to the next point*/ + griddata[i*ncols+j]=data_value; + } + } + } + } + if (debug) printf("\b\b\b\b\b\b\b%5.2lf %%\n",100.0); + + /*Assign output pointers:*/ + *pgriddata=griddata; + *px_m=x_m; + *py_m=y_m; +} Index: /issm/trunk/src/c/modules/InterpFromMeshToGridx/InterpFromMeshToGridx.h =================================================================== --- /issm/trunk/src/c/modules/InterpFromMeshToGridx/InterpFromMeshToGridx.h (revision 3904) +++ /issm/trunk/src/c/modules/InterpFromMeshToGridx/InterpFromMeshToGridx.h (revision 3904) @@ -0,0 +1,12 @@ +/*!\file InterpFromMeshToGridx.h + * \brief: header file for Data interpolation routines. + */ + +#ifndef _INTERPFROMMESHTOGRIDX_H +#define _INTERPFROMMESHTOGRIDX_H + +#include "../toolkits/toolkits.h" + +void InterpFromMeshToGridx(double** px_m,double** py_m,double** pgriddata,double* index_mesh, double* x_mesh, double* y_mesh, int nods,int nels, double* data_mesh, int data_length, double cornereast,double cornernorth,double xposting,double yposting,int nlines,int ncols,double default_value); + +#endif /* _INTERPFROMMESHTOGRIDX_H */ Index: /issm/trunk/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp =================================================================== --- /issm/trunk/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp (revision 3904) +++ /issm/trunk/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp (revision 3904) @@ -0,0 +1,118 @@ +/*!\file InterpFromMeshToMesh2dx + */ + +#include "./InterpFromMeshToMesh2dx.h" + +#include "../shared/shared.h" +#include "../include/include.h" +#include "../toolkits/toolkits.h" + +#include "../Bamgx/objects/BamgObjects.h" + +using namespace bamg; +using namespace std; + +int InterpFromMeshToMesh2dx(double** pdata_interp,double* index_data,double* x_data,double* y_data,int nods_data,int nels_data, + double* data,int data_rows,int data_cols,double* x_interp,double* y_interp,int nods_interp,double default_value){ + + /*Output*/ + double* data_interp=NULL; + int noerr; + + /*Intermediary*/ + R2 r; + I2 I; + int i,j,k; + int it; + int i0,i1,i2; + double areacoord[3]; + double aa,bb; + double data_value; + Icoor2 dete[3]; + int verbose=0; + + /*Checks*/ + if (data_cols<=0){ + ISSMERROR("data provided has a negative number of columns"); + } + if (data_rows!=nods_data && data_rows!=nels_data){ + ISSMERROR("data provided should have either %i or %i lines (not %i)",nods_data,nels_data,data_rows); + } + + /*Initialize output*/ + if (verbose) printf("Initializing output vector\n"); + data_interp=(double*)xmalloc(nods_interp*data_cols*sizeof(double)); + + // read background mesh + if (verbose) printf("Reading mesh\n"); + Triangles Th(index_data,x_data,y_data,nods_data,nels_data); + Th.ReMakeTriangleContainingTheVertex(); + + //Loop over output nodes + if (verbose) printf("Loop over the nodes\n"); + for(i=0;i0){ + //Area coordinate + areacoord[0]= (double) dete[0]/ tb.det; + areacoord[1]= (double) dete[1] / tb.det; + areacoord[2]= (double) dete[2] / tb.det; + //3 vertices of the triangle + i0=Th.Number(tb[0]); + i1=Th.Number(tb[1]); + i2=Th.Number(tb[2]); + //triangle number + it=Th.Number(tb); + } + + //external point + else { + //Get closest adjacent triangle (inside the mesh) + TriangleAdjacent ta=CloseBoundaryEdge(I,&tb,aa,bb).Adj(); + int k=ta; + Triangle &tc=*(Triangle*)ta; + //Area coordinate + areacoord[VerticesOfTriangularEdge[k][1]] = aa; + areacoord[VerticesOfTriangularEdge[k][0]] = bb; + areacoord[OppositeVertex[k]] = 1 - aa -bb; + //3 vertices of the triangle + i0=Th.Number(tc[0]); + i1=Th.Number(tc[1]); + i2=Th.Number(tc[2]); + //triangle number + it=Th.Number(tc); + } + + /*Last step, P1 interpolation*/ + if (data_rows==nods_data){ + for (j=0;j=nels_data){ + ISSMERROR("Triangle number %i not in [0 %i], because not correctly implemented yet... interpolate on grid first",it,nels_data); + } + data_interp[i*data_cols+j]=data[data_cols*it+j]; + } + } + + } + + /*Assign output pointers:*/ + if (verbose) printf("Assigning output\n"); + *pdata_interp=data_interp; + + /*No error return*/ + return noerr; + +} Index: /issm/trunk/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.h =================================================================== --- /issm/trunk/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.h (revision 3904) +++ /issm/trunk/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.h (revision 3904) @@ -0,0 +1,14 @@ +/*!\file: InterpFromMeshToMesh2dx.h + * \brief header file for Bamg module + */ + +#ifndef _INTERPFROMMESHTOMESH2DX_H +#define _INTERPFROMMESHTOMESH2DX_H + +#include "../objects/objects.h" + +/* local prototypes: */ +int InterpFromMeshToMesh2dx(double** pdata_interp,double* index_data,double* x_data,double* y_data,int nods_data,int nels_data, + double* data,int data_rows,int data_cols,double* x_interp,double* y_interp,int nods_interp,double default_value); + +#endif Index: /issm/trunk/src/c/modules/InterpFromMeshToMesh3dx/InterpFromMeshToMesh3dx.cpp =================================================================== --- /issm/trunk/src/c/modules/InterpFromMeshToMesh3dx/InterpFromMeshToMesh3dx.cpp (revision 3904) +++ /issm/trunk/src/c/modules/InterpFromMeshToMesh3dx/InterpFromMeshToMesh3dx.cpp (revision 3904) @@ -0,0 +1,138 @@ +/*!\file: InterpFromMeshToMesh3dx.cpp + * \brief "c" core code for interpolating values from a structured grid. + */ + +#include "./InterpFromMeshToMesh3dx.h" +#include "../shared/shared.h" +#include "../include/include.h" + +int InterpFromMeshToMesh3dx( Vec* pdata_prime,double* index_data, double* x_data, double* y_data, double* z_data, int nods_data,int nels_data, double* data, int data_length, double* x_prime, double* y_prime, double* z_prime, int nods_prime,double default_value) { + + /*Output*/ + Vec data_prime=NULL; + + /*Intermediary*/ + int i,j; + int interpolation_type; + bool debug; + double area; + double area_1,area_2,area_3; + double zeta,bed,surface; + double data_value; + double x_prime_min,x_prime_max; + double y_prime_min,y_prime_max; + double x_tria_min,y_tria_min; + double x_tria_max,y_tria_max; + + /*some checks*/ + if (nels_data<1 || nods_data<6 || nods_prime==0){ + ISSMERROR("nothing to be done according to the mesh given in input"); + } + + /*Set debug to 1 if there are lots of elements*/ + debug=(bool)((double)nels_data*(double)nods_prime >= pow((double)10,(double)9)); + + /*figure out what kind of interpolation is needed*/ + if (data_length==nods_data){ + interpolation_type=1; + } + else if (data_length==nels_data){ + interpolation_type=2; + } + else{ + ISSMERROR("length of vector data not supported yet. It should be of length (number of nodes) or (number of elements)!"); + } + + /*Get prime mesh extrema coordinates*/ + x_prime_min=x_prime[0]; x_prime_max=x_prime[0];y_prime_min=y_prime[0]; y_prime_max=y_prime[0]; + for (i=1;ix_prime_max) x_prime_max=x_prime[i]; + if (y_prime[i]y_prime_max) y_prime_max=y_prime[i]; + } + + /*Initialize output*/ + data_prime=NewVec(nods_prime); + for (i=0;ix_tria_max) x_tria_max=x_data[(int)index_data[6*i+j]-1]; + if(y_data[(int)index_data[6*i+j]-1]y_tria_max) y_tria_max=y_data[(int)index_data[6*i+j]-1]; + } + + /*if there is no point inside the domain, go to next iteration*/ + if ( x_prime_max < x_tria_min ) continue; + if ( x_prime_min > x_tria_max ) continue; + if ( y_prime_max < y_tria_min ) continue; + if ( y_prime_min > y_tria_max ) continue; + + /*get area of the current element (Jacobian = 2 * area)*/ + //area =x2 * y3 - y2*x3 + x1 * y2 - y1 * x2 + x3 * y1 - y3 * x1; + area=x_data[(int)index_data[6*i+1]-1]*y_data[(int)index_data[6*i+2]-1]-y_data[(int)index_data[6*i+1]-1]*x_data[(int)index_data[6*i+2]-1] + + x_data[(int)index_data[6*i+0]-1]*y_data[(int)index_data[6*i+1]-1]-y_data[(int)index_data[6*i+0]-1]*x_data[(int)index_data[6*i+1]-1] + + x_data[(int)index_data[6*i+2]-1]*y_data[(int)index_data[6*i+0]-1]-y_data[(int)index_data[6*i+2]-1]*x_data[(int)index_data[6*i+0]-1]; + + /*loop over the prime nodes*/ + for (j=0;j x_tria_max ) continue; + if ( y_prime[j] < y_tria_min ) continue; + if ( y_prime[j] > y_tria_max ) continue; + + /*Get first area coordinate = det(x-x3 x2-x3 ; y-y3 y2-y3)/area*/ + area_1=((x_prime[j]-x_data[(int)index_data[6*i+2]-1])*(y_data[(int)index_data[6*i+1]-1]-y_data[(int)index_data[6*i+2]-1]) + - (y_prime[j]-y_data[(int)index_data[6*i+2]-1])*(x_data[(int)index_data[6*i+1]-1]-x_data[(int)index_data[6*i+2]-1]))/area; + /*Get second area coordinate =det(x1-x3 x-x3 ; y1-y3 y-y3)/area*/ + area_2=((x_data[(int)index_data[6*i+0]-1]-x_data[(int)index_data[6*i+2]-1])*(y_prime[j]-y_data[(int)index_data[6*i+2]-1]) + - (y_data[(int)index_data[6*i+0]-1]-y_data[(int)index_data[6*i+2]-1])*(x_prime[j]-x_data[(int)index_data[6*i+2]-1]))/area; + /*Get third area coordinate = 1-area1-area2*/ + area_3=1-area_1-area_2; + + /*is the current point in the current 2d element?*/ + if (area_1>=0 && area_2>=0 && area_3>=0){ + + /*compute bottom and top height of the element at this 2d position*/ + bed =area_1*z_data[(int)index_data[6*i+0]-1]+area_2*z_data[(int)index_data[6*i+1]-1]+area_3*z_data[(int)index_data[6*i+2]-1]; + surface=area_1*z_data[(int)index_data[6*i+3]-1]+area_2*z_data[(int)index_data[6*i+4]-1]+area_3*z_data[(int)index_data[6*i+5]-1]; + + /*Compute zeta*/ + zeta=2*(z_prime[j]-bed)/(surface-bed)-1; + + if (zeta >=-1 && zeta<=1){ + if (interpolation_type==1){ + /*nodal interpolation*/ + data_value=(1-zeta)/2*(area_1*data[(int)index_data[6*i+0]-1]+area_2*data[(int)index_data[6*i+1]-1]+area_3*data[(int)index_data[6*i+2]-1]) + + (1+zeta)/2*(area_1*data[(int)index_data[6*i+3]-1]+area_2*data[(int)index_data[6*i+4]-1]+area_3*data[(int)index_data[6*i+5]-1]); + } + else{ + /*element interpolation*/ + data_value=data[i]; + } + if (isnan(data_value)) data_value=default_value; + + /*insert value and go to the next point*/ + VecSetValue(data_prime,j,data_value,INSERT_VALUES); + } + } + } + } + if (debug) printf("\b\b\b\b\b\b\b%5.2lf %%\n",100.0); + + /*Assign output pointers:*/ + *pdata_prime=data_prime; +} Index: /issm/trunk/src/c/modules/InterpFromMeshToMesh3dx/InterpFromMeshToMesh3dx.h =================================================================== --- /issm/trunk/src/c/modules/InterpFromMeshToMesh3dx/InterpFromMeshToMesh3dx.h (revision 3904) +++ /issm/trunk/src/c/modules/InterpFromMeshToMesh3dx/InterpFromMeshToMesh3dx.h (revision 3904) @@ -0,0 +1,13 @@ +/*!\file InterpFromMeshToMesh3dx.h + * \brief: header file for Data interpolation routines. + */ + +#ifndef _INTERPFROMMESHTOMESH3DX_H +#define _INTERPFROMMESHTOMESH3DX_H + +#include "../toolkits/toolkits.h" + +int InterpFromMeshToMesh3dx( Vec* pdata_prime,double* index_data, double* x_data, double* y_data, double* z_data, int nods_data,int nels_data, double* data, int data_length, double* x_prime, double* y_prime, double* z_prime, int nods_prime,double default_value); + +#endif /* _INTERPFROMMESHTOMESH3DX_H */ + Index: /issm/trunk/src/c/modules/MassFluxx/MassFluxx.cpp =================================================================== --- /issm/trunk/src/c/modules/MassFluxx/MassFluxx.cpp (revision 3904) +++ /issm/trunk/src/c/modules/MassFluxx/MassFluxx.cpp (revision 3904) @@ -0,0 +1,53 @@ +/*!\file MassFluxx + * \brief: compute mass flux along a profile. + */ + +#include "./MassFluxx.h" + +#include "../shared/shared.h" +#include "../include/include.h" +#include "../toolkits/toolkits.h" +#include "../EnumDefinitions/EnumDefinitions.h" + +void MassFluxx(double* pmass_flux, DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads,DataSet* materials,Parameters* parameters, double* segments,int num_segments,double* ug){ + + int i,j; + extern int num_procs; + extern int my_rank; + + Element* element=NULL; + int element_id; + + /*output: */ + double mass_flux=0; + double all_mass_flux=0; + + /*First, get elements and loads configured: */ + elements->Configure(elements,loads, nodes,vertices, materials,parameters); + nodes->Configure(elements, loads, nodes,vertices, materials,parameters); + loads->Configure(elements, loads, nodes,vertices, materials,parameters); + parameters->Configure(elements,loads, nodes,vertices, materials,parameters); + + /*Go through segments, and then elements, and figure out which elements belong to a segment. + * When we find one, use the element to compute the mass flux on the segment: */ + + for(i=0;iSize();j++){ + element=(Element*)elements->GetObjectByOffset(j); + if (element->Id()==element_id){ + /*We found the element which owns this segment, use it to compute the mass flux: */ + mass_flux+=element->MassFlux(segments+5*i+0,ug); + break; + } + } + } + + #ifdef _PARALLEL_ + MPI_Allreduce ( (void*)&mass_flux,(void*)&all_mass_flux,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); + mass_flux=all_mass_flux; + #endif + + /*Assign output pointers: */ + *pmass_flux=mass_flux; +} Index: /issm/trunk/src/c/modules/MassFluxx/MassFluxx.h =================================================================== --- /issm/trunk/src/c/modules/MassFluxx/MassFluxx.h (revision 3904) +++ /issm/trunk/src/c/modules/MassFluxx/MassFluxx.h (revision 3904) @@ -0,0 +1,16 @@ +/*!\file: MassFluxx.h + * \brief header file for degree of freedoms distribution routines. + */ + +#ifndef _MASSFLUXX_H +#define _MASSFLUXX_H + +#include "../DataSet/DataSet.h" +#include "../objects/objects.h" + +/* local prototypes: */ +void MassFluxx(double* pmass_flux, DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads,DataSet* materials, Parameters* parameters,double* segments,int num_segments,double* ug); + + +#endif /* _MASSFLUXX_H */ + Index: /issm/trunk/src/c/modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp =================================================================== --- /issm/trunk/src/c/modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp (revision 3904) +++ /issm/trunk/src/c/modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp (revision 3904) @@ -0,0 +1,48 @@ +/*!\file Mergesolutionfromftogx + * \brief merge solution back from f set into g set + */ + +#include "./Mergesolutionfromftogx.h" + +void Mergesolutionfromftogx( Vec* pug, Vec uf, Mat Gmn, Vec ys, NodeSets* nodesets){ + + /*output: */ + Vec ug=NULL; + int ug_local_size; + + /*intermediary*/ + Vec u_m=NULL; + Vec u_n=NULL; + int u_n_local_size; + int Gmn_localm,Gmn_localn; + + + /*Merge f set back into n set: */ + u_n=NewVec(nodesets->GetNSize()); + VecMerge(u_n,uf,nodesets->GetPV_F(),nodesets->GetFSize()); + + /*merge single point constraints into n set: */ + if(ys){ + VecMerge(u_n,ys,nodesets->GetPV_S(),nodesets->GetSSize()); + } + + /* Merge n_set and m_set to g_set:*/ + ug=NewVec(nodesets->GetGSize()); + VecMerge(ug,u_n,nodesets->GetPV_N(),nodesets->GetNSize()); + + /*Merge m set into g set: */ + if(nodesets->GetMSize()){ + MatGetLocalSize(Gmn,&Gmn_localm,&Gmn_localn); + u_m=NewVecFromLocalSize(Gmn_localm); + MatMultPatch(Gmn,u_n,u_m); + VecMerge(ug,u_m,nodesets->GetPV_M(),nodesets->GetMSize()); + } + + + /*Free ressources:*/ + VecFree(&u_m); + VecFree(&u_n); + + /*Assign correct pointer*/ + *pug=ug; +} Index: /issm/trunk/src/c/modules/Mergesolutionfromftogx/Mergesolutionfromftogx.h =================================================================== --- /issm/trunk/src/c/modules/Mergesolutionfromftogx/Mergesolutionfromftogx.h (revision 3904) +++ /issm/trunk/src/c/modules/Mergesolutionfromftogx/Mergesolutionfromftogx.h (revision 3904) @@ -0,0 +1,14 @@ +/*!\file: Mergesolutionfromftogx.h + * \brief merge solution back from f set into g set + */ + +#ifndef _MERGESOLUTIONFROMFTOGX_H +#define _MERGESOLUTIONFROMFTOGX_H + +#include "../objects/objects.h" + +/* local prototypes: */ +void Mergesolutionfromftogx( Vec* pug, Vec uf, Mat Gmn, Vec ys, NodeSets* nodesets); + +#endif /* _MERGESOLUTIONFROMFTOGX_H */ + Index: /issm/trunk/src/c/modules/MeshPartitionx/MeshPartitionx.cpp =================================================================== --- /issm/trunk/src/c/modules/MeshPartitionx/MeshPartitionx.cpp (revision 3904) +++ /issm/trunk/src/c/modules/MeshPartitionx/MeshPartitionx.cpp (revision 3904) @@ -0,0 +1,113 @@ +/*!\file: MeshPartition.cpp + * \brief partition elements and grids across a cluster of size numprocs. + */ + +#include "./MeshPartitionx.h" +#include "../shared/shared.h" +#include "../include/include.h" +#include "../toolkits/toolkits.h" +#include "../EnumDefinitions/EnumDefinitions.h" + +int MeshPartitionx(int** pepart, int** pnpart, int numberofelements,int numberofgrids,double* elements, + int numberofelements2d,int numberofgrids2d,double* elements2d,int numlayers,int elements_width, char* meshtype,int num_procs){ + + int noerr=1; + int i,j; + + /*Metis partitioning: */ + int* epart=NULL; + int* npart=NULL; + int* index=NULL; + + int* epart2d=NULL; + int* npart2d=NULL; + int* index2d=NULL; + int count=0; + + int etype=1; //tria mesh see metis/Programs/Io.c + int etype2d=1; //tria mesh see metis/Programs/Io.c + int numflag=0; + int edgecut=1; + + if(strcmp(meshtype,"2d")==0){ + epart=(int*)xmalloc(numberofelements*sizeof(int)); + npart=(int*)xmalloc(numberofgrids*sizeof(int)); + index=(int*)xmalloc(elements_width*numberofelements*sizeof(int)); + for (i=0;i1){ + METIS_PartMeshNodal(&numberofelements,&numberofgrids, index, &etype, &numflag, &num_procs, &edgecut, epart, npart); + } + else if (num_procs==1){ + /*METIS does not know how to deal with one cpu only!*/ + for (i=0;i1){ + METIS_PartMeshNodal(&numberofelements2d,&numberofgrids2d, index2d, &etype2d, &numflag, &num_procs, &edgecut, epart2d, npart2d); + } + else if (num_procs==1){ + /*METIS does not know how to deal with one cpu only!*/ + for (i=0;i