Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 7639) +++ /issm/trunk/src/c/Makefile.am (revision 7640) @@ -395,4 +395,8 @@ ./modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp\ ./modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp\ + ./modules/ModelProcessorx/Hydrology/UpdateElementsHydrology.cpp\ + ./modules/ModelProcessorx/Hydrology/CreateNodesHydrology.cpp\ + ./modules/ModelProcessorx/Hydrology/CreateConstraintsHydrology.cpp\ + ./modules/ModelProcessorx/Hydrology/CreateLoadsHydrology.cpp\ ./modules/ModelProcessorx/Melting/UpdateElementsMelting.cpp\ ./modules/ModelProcessorx/Melting/CreateNodesMelting.cpp\ @@ -420,4 +424,7 @@ ./modules/OutputResultsx/OutputResultsx.h\ ./modules/OutputResultsx/OutputResultsx.cpp\ + ./modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.h\ + ./modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.cpp\ + ./modules/PointCloudFindNeighborsx/PointCloudFindNeighborsxt.cpp\ ./modules/MinVelx/MinVelx.h\ ./modules/MinVelx/MinVelx.cpp\ @@ -970,4 +977,8 @@ ./modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp\ ./modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp\ + ./modules/ModelProcessorx/Hydrology/UpdateElementsHydrology.cpp\ + ./modules/ModelProcessorx/Hydrology/CreateNodesHydrology.cpp\ + ./modules/ModelProcessorx/Hydrology/CreateConstraintsHydrology.cpp\ + ./modules/ModelProcessorx/Hydrology/CreateLoadsHydrology.cpp\ ./modules/ModelProcessorx/Melting/UpdateElementsMelting.cpp\ ./modules/ModelProcessorx/Melting/CreateNodesMelting.cpp\ @@ -1171,4 +1182,6 @@ ./solutions/surfaceslope_core.cpp\ ./solutions/bedslope_core.cpp\ + ./solutions/hydrology_core.cpp\ + ./solutions/hydrology_core_step.cpp\ ./solutions/transient2d_core.cpp\ ./solutions/transient3d_core.cpp\ @@ -1182,5 +1195,5 @@ ./solvers/solver_linear.cpp\ ./solvers/solver_adjoint_linear.cpp\ - ./solvers/solver_diagnostic_nonlinear.cpp\ + ./solvers/solver_nonlinear.cpp\ ./solvers/solver_stokescoupling_nonlinear.cpp\ ./solvers/solver_thermal_nonlinear.cpp\ Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp (revision 7639) +++ /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp (revision 7640) @@ -72,4 +72,11 @@ UpdateElementsThermal(elements,iomodel,iomodel_handle,analysis_counter,analysis_type); break; + + case HydrologyAnalysisEnum: + CreateNodesHydrology(pnodes, iomodel,iomodel_handle); + CreateConstraintsHydrology(pconstraints,iomodel,iomodel_handle); + CreateLoadsHydrology(ploads,iomodel,iomodel_handle); + UpdateElementsHydrology(elements,iomodel,iomodel_handle,analysis_counter,analysis_type); + break; case MeltingAnalysisEnum: Index: /issm/trunk/src/c/modules/ModelProcessorx/DistributeNumDofs.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/DistributeNumDofs.cpp (revision 7639) +++ /issm/trunk/src/c/modules/ModelProcessorx/DistributeNumDofs.cpp (revision 7640) @@ -75,4 +75,7 @@ numdofs=1; } + else if (analysis_type==HydrologyAnalysisEnum){ + numdofs=1; + } else if (analysis_type==MeltingAnalysisEnum){ numdofs=1; Index: /issm/trunk/src/c/modules/ModelProcessorx/Hydrology/CreateConstraintsHydrology.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Hydrology/CreateConstraintsHydrology.cpp (revision 7640) +++ /issm/trunk/src/c/modules/ModelProcessorx/Hydrology/CreateConstraintsHydrology.cpp (revision 7640) @@ -0,0 +1,54 @@ +/* + * CreateConstraintsHydrology.c: + */ + +#include "../../../Container/Container.h" +#include "../../../toolkits/toolkits.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../io/io.h" +#include "../../../objects/objects.h" +#include "../../../shared/shared.h" +#include "../ModelProcessorx.h" + +void CreateConstraintsHydrology(Constraints** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){ + + /*Intermediary*/ + int i; + int count; + + /*Output*/ + Constraints *constraints = NULL; + Spc *spc = NULL; + + /*Recover pointer: */ + constraints=*pconstraints; + + /*Create constraints if they do not exist yet*/ + if(!constraints) constraints = new Constraints(ConstraintsEnum); + + /*Fetch data: */ + IoModelFetchData(&iomodel->spcwatercolumn,NULL,NULL,iomodel_handle,"spcwatercolumn"); + + /*Initialize counter*/ + count=0; + + /*Create spcs from x,y,z, as well as the spc values on those spcs: */ + for (i=0;inumberofvertices;i++){ + + /*keep only this partition's nodes:*/ + if((iomodel->my_vertices[i])){ + + if ((int)iomodel->spcwatercolumn[2*i]){ + + constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,*(iomodel->spcwatercolumn+2*i+1),HydrologyAnalysisEnum)); + count++; + } + } //if((my_vertices[i])) + } + + /*Free data: */ + xfree((void**)&iomodel->spcwatercolumn); + + /*Assign output pointer: */ + *pconstraints=constraints; +} Index: /issm/trunk/src/c/modules/ModelProcessorx/Hydrology/CreateLoadsHydrology.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Hydrology/CreateLoadsHydrology.cpp (revision 7640) +++ /issm/trunk/src/c/modules/ModelProcessorx/Hydrology/CreateLoadsHydrology.cpp (revision 7640) @@ -0,0 +1,29 @@ +/*! \file CreateLoadsHydrology.c: + */ + +#include "../../../Container/Container.h" +#include "../../../toolkits/toolkits.h" +#include "../../../io/io.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../objects/objects.h" +#include "../../../shared/shared.h" +#include "../../../include/include.h" +#include "../ModelProcessorx.h" + +void CreateLoadsHydrology(Loads** ploads, IoModel* iomodel,ConstDataHandle iomodel_handle){ + + /*Intermediary*/ + int i; + + /*DataSet*/ + Loads* loads = NULL; + + /*Recover pointer: */ + loads=*ploads; + + /*Create loads if they do not exist yet*/ + if(!loads) loads = new Loads(LoadsEnum); + + /*Assign output pointer: */ + *ploads=loads; +} Index: /issm/trunk/src/c/modules/ModelProcessorx/Hydrology/CreateNodesHydrology.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Hydrology/CreateNodesHydrology.cpp (revision 7640) +++ /issm/trunk/src/c/modules/ModelProcessorx/Hydrology/CreateNodesHydrology.cpp (revision 7640) @@ -0,0 +1,61 @@ +/* + * CreateNodesHydrology.c: + */ + +#include "../../../Container/Container.h" +#include "../../../toolkits/toolkits.h" +#include "../../../io/io.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../objects/objects.h" +#include "../../../include/include.h" +#include "../../../shared/shared.h" +#include "../../MeshPartitionx/MeshPartitionx.h" +#include "../ModelProcessorx.h" + +void CreateNodesHydrology(Nodes** pnodes, IoModel* iomodel,ConstDataHandle iomodel_handle){ + + /*Intermediary*/ + int i; + bool continuous_galerkin=true; + + /*DataSets: */ + Nodes* nodes = NULL; + + /*Recover pointer: */ + nodes=*pnodes; + + /*Create nodes if they do not exist yet*/ + if(!nodes) nodes = new Nodes(NodesEnum); + + /*Continuous Galerkin partition of nodes: */ + NodesPartitioning(&iomodel->my_nodes,iomodel->my_elements,iomodel->my_vertices,iomodel,iomodel_handle,continuous_galerkin); + + /*Create nodes and vertices: */ + IoModelFetchData(&iomodel->gridonbed,NULL,NULL,iomodel_handle,"gridonbed"); + IoModelFetchData(&iomodel->gridonsurface,NULL,NULL,iomodel_handle,"gridonsurface"); + IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet"); + IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf"); + IoModelFetchData(&iomodel->vertices_type,NULL,NULL,iomodel_handle,"vertices_type"); + IoModelFetchData(&iomodel->gridonwater,NULL,NULL,iomodel_handle,"gridonwater"); + + for (i=0;inumberofvertices;i++){ + + if(iomodel->my_vertices[i]){ + + /*Add node to nodes dataset: */ + nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i+1,i,iomodel,HydrologyAnalysisEnum)); + + } + } + + /*Clean fetched data: */ + xfree((void**)&iomodel->gridonbed); + xfree((void**)&iomodel->gridonsurface); + xfree((void**)&iomodel->gridonicesheet); + xfree((void**)&iomodel->gridonwater); + xfree((void**)&iomodel->gridoniceshelf); + xfree((void**)&iomodel->vertices_type); + + /*Assign output pointer: */ + *pnodes=nodes; +} Index: /issm/trunk/src/c/modules/ModelProcessorx/Hydrology/UpdateElementsHydrology.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Hydrology/UpdateElementsHydrology.cpp (revision 7640) +++ /issm/trunk/src/c/modules/ModelProcessorx/Hydrology/UpdateElementsHydrology.cpp (revision 7640) @@ -0,0 +1,81 @@ +/* + * UpdateElementsHydrology: + */ + +#include "../../../Container/Container.h" +#include "../../../toolkits/toolkits.h" +#include "../../../io/io.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../objects/objects.h" +#include "../../../shared/shared.h" +#include "../../MeshPartitionx/MeshPartitionx.h" +#include "../../../include/include.h" +#include "../ModelProcessorx.h" + +void UpdateElementsHydrology(Elements* elements, IoModel* iomodel,ConstDataHandle iomodel_handle,int analysis_counter,int analysis_type){ + + /*Intermediary*/ + int i; + int counter; + Element* element=NULL; + + /*Fetch data needed: */ + IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); + IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness"); + IoModelFetchData(&iomodel->surface,NULL,NULL,iomodel_handle,"surface"); + IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed"); + IoModelFetchData(&iomodel->drag_coefficient,NULL,NULL,iomodel_handle,"drag_coefficient"); + IoModelFetchData(&iomodel->drag_p,NULL,NULL,iomodel_handle,"drag_p"); + IoModelFetchData(&iomodel->drag_q,NULL,NULL,iomodel_handle,"drag_q"); + IoModelFetchData(&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf"); + IoModelFetchData(&iomodel->elementonbed,NULL,NULL,iomodel_handle,"elementonbed"); + IoModelFetchData(&iomodel->elementonsurface,NULL,NULL,iomodel_handle,"elementonsurface"); + IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater"); + IoModelFetchData(&iomodel->elements_type,NULL,NULL,iomodel_handle,"elements_type"); + IoModelFetchData(&iomodel->rheology_B,NULL,NULL,iomodel_handle,"rheology_B"); + IoModelFetchData(&iomodel->rheology_n,NULL,NULL,iomodel_handle,"rheology_n"); + IoModelFetchData(&iomodel->pressure,NULL,NULL,iomodel_handle,"pressure"); + IoModelFetchData(&iomodel->temperature,NULL,NULL,iomodel_handle,"temperature"); + IoModelFetchData(&iomodel->melting_rate,NULL,NULL,iomodel_handle,"melting_rate"); + IoModelFetchData(&iomodel->watercolumn,NULL,NULL,iomodel_handle,"watercolumn"); + IoModelFetchData(&iomodel->geothermalflux,NULL,NULL,iomodel_handle,"geothermalflux"); + IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx"); + IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy"); + IoModelFetchData(&iomodel->vz,NULL,NULL,iomodel_handle,"vz"); + + /*Update elements: */ + counter=0; + for (i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type); //we need i to index into elements. + counter++; + } + } + + cleanup_and_return: + + /*Free data: */ + xfree((void**)&iomodel->elements); + xfree((void**)&iomodel->thickness); + xfree((void**)&iomodel->surface); + xfree((void**)&iomodel->bed); + xfree((void**)&iomodel->drag_coefficient); + xfree((void**)&iomodel->drag_p); + xfree((void**)&iomodel->drag_q); + xfree((void**)&iomodel->elementoniceshelf); + xfree((void**)&iomodel->elementonbed); + xfree((void**)&iomodel->elementonsurface); + xfree((void**)&iomodel->elementonwater); + xfree((void**)&iomodel->elements_type); + xfree((void**)&iomodel->rheology_B); + xfree((void**)&iomodel->rheology_n); + xfree((void**)&iomodel->pressure); + xfree((void**)&iomodel->temperature); + xfree((void**)&iomodel->melting_rate); + xfree((void**)&iomodel->watercolumn); + xfree((void**)&iomodel->geothermalflux); + xfree((void**)&iomodel->vx); + xfree((void**)&iomodel->vy); + xfree((void**)&iomodel->vz); +} Index: /issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.h =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 7639) +++ /issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 7640) @@ -61,4 +61,11 @@ void UpdateElementsThermal(Elements* elements,IoModel* iomodel,ConstDataHandle iomodel_handle,int analysis_counter,int analysis_type); +/*hydrology:*/ +void CreateNodesHydrology(Nodes** pnodes,IoModel* iomodel,ConstDataHandle iomodel_handle); +void CreateConstraintsHydrology(Constraints** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); +void CreateLoadsHydrology(Loads** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); +void UpdateElementsHydrology(Elements* elements,IoModel* iomodel,ConstDataHandle iomodel_handle,int analysis_counter,int analysis_type); + + /*melting:*/ void CreateNodesMelting(Nodes** pnodes,IoModel* iomodel,ConstDataHandle iomodel_handle); Index: /issm/trunk/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.cpp =================================================================== --- /issm/trunk/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.cpp (revision 7640) +++ /issm/trunk/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.cpp (revision 7640) @@ -0,0 +1,38 @@ +/*! \file PointCloudFindNeighborsx.c + */ + +#include "./PointCloudFindNeighborsx.h" + +int PointCloudFindNeighborsx( Vec* pflags,double* x, double* y, int nods, double mindistance,double multithread){ + + /*output: */ + Vec flags=NULL; + flags=NewVec(nods); + + /*threading: */ + PointCloudFindNeighborsThreadStruct gate; + int num=1; + + #ifdef _MULTITHREADING_ + num=_NUMTHREADS_; + #endif + + if(!multithread)num=1; + + /*initialize thread parameters: */ + gate.x=x; + gate.y=y; + gate.nods=nods; + gate.mindistance=mindistance; + gate.flags=flags; + + /*launch the thread manager with InterpFromGridToMeshxt as a core: */ + LaunchThread(PointCloudFindNeighborsxt,(void*)&gate,num); + + /*Assemble vector: */ + VecAssemblyBegin(flags); + VecAssemblyEnd(flags); + + /*Assign output pointers: */ + *pflags=flags; +} Index: /issm/trunk/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.h =================================================================== --- /issm/trunk/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.h (revision 7640) +++ /issm/trunk/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.h (revision 7640) @@ -0,0 +1,31 @@ +/* + PointCloudFindNeighborsx.h +*/ + + +#ifndef _POINTCLOUDFLAGNEIGHBORSX_H +#define _POINTCLOUDFLAGNEIGHBORSX_H + +#include "../../shared/shared.h" +#include "../../objects/objects.h" + +/* local prototypes: */ +int PointCloudFindNeighborsx( Vec* pflags,double* x, double* y, int nods, double mindistance,double multithread); + +/*threading: */ +typedef struct{ + + double* x; + double* y; + int nods; + double mindistance; + Vec flags; + + +} PointCloudFindNeighborsThreadStruct; + +void* PointCloudFindNeighborsxt(void* vPointCloudFindNeighborsThreadStruct); + + +#endif /* _POINTCLOUDFLAGNEIGHBORSX_H */ + Index: /issm/trunk/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsxt.cpp =================================================================== --- /issm/trunk/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsxt.cpp (revision 7640) +++ /issm/trunk/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsxt.cpp (revision 7640) @@ -0,0 +1,78 @@ +/*!\file: PointCloudFindNeighborst.cpp + * \brief thread core for PointCloudFindNeighborst code + */ + +#include "./PointCloudFindNeighborsx.h" +#include "../../shared/shared.h" + +void* PointCloudFindNeighborsxt(void* vpthread_handle){ + + /*gate variables :*/ + PointCloudFindNeighborsThreadStruct* gate=NULL; + pthread_handle* handle=NULL; + int my_thread; + int num_threads; + + double* x; + double* y; + int nods; + double mindistance; + Vec flags; + + /*recover handle and gate: */ + handle=(pthread_handle*)vpthread_handle; + gate=(PointCloudFindNeighborsThreadStruct*)handle->gate; + my_thread=handle->id; + num_threads=handle->num; + + /*recover parameters :*/ + x=gate->x; + y=gate->y; + nods=gate->nods; + mindistance=gate->mindistance; + flags=gate->flags; + + /*intermediary: */ + int i,j; + int i0,i1; + double distance; + bool* already=NULL; + + /*allocate: */ + already=(bool*)xcalloc(nods,sizeof(bool)); + + /*partition loop across threads: */ + PartitionRange(&i0,&i1,nods,num_threads,my_thread); + + /*Loop over the nodes*/ + if (my_thread==0) printf(" interpolation progress: %5.2lf %%",0.0); + + for (i=i0;i error*/ if( v && Norme1(v->r - vertices[i].r) < eps ){ - _error_("two points of the geometry are very closed to each other"); + //_error_("two points of the geometry are very closed to each other"); + printf("%g %g\n",v->ReferenceNumber,vertices[i].ReferenceNumber); } Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 7639) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 7640) @@ -405,4 +405,7 @@ Ke=CreateKMatrixPrognostic(); break; + case HydrologyAnalysisEnum: + Ke=CreateKMatrixHydrology(); + break; case BalancedthicknessAnalysisEnum: Ke=CreateKMatrixBalancedthickness(); @@ -1182,4 +1185,85 @@ } /*}}}*/ +/*FUNCTION Tria::CreateKMatrixHydrology {{{1*/ +ElementMatrix* Tria::CreateKMatrixHydrology(void){ + + /*Constants*/ + const int numdof=NDOF1*NUMVERTICES; + + /*Intermediaries */ + int i,j,ig; + double Jdettria,dt; + double xyz_list[NUMVERTICES][3]; + double L[NUMVERTICES]; + double dL[2][NUMVERTICES]; + double dLx[NUMVERTICES]; + double dLy[NUMVERTICES]; + double K[2]; + double Ke_gg_gaussian1[numdof][numdof] ={0.0}; + double Ke_gg_gaussian2[numdof][numdof] ={0.0}; + double Ke_gg_gaussian3[numdof][numdof] ={0.0}; + double gamma; + double DL_scalar1; + double DL_scalar2; + double DL_scalar3; + GaussTria *gauss=NULL; + + /*Initialize Element matrix*/ + ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + this->parameters->FindParam(&dt,DtEnum); + + /*retrieve material parameters: */ + gamma=matpar->GetGamma(); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(2); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); + GetNodalFunctions(&L[0],gauss); + GetNodalFunctionsDerivatives(&dL[0][0],&xyz_list[0][0],gauss); + + for(i=0;iweight*Jdettria; + TripleMultiply( &L[0],1,numdof,1, + &DL_scalar1,1,1,0, + &L[0],1,numdof,0, + &Ke_gg_gaussian1[0][0],0); + } + + if(dt){ + DL_scalar2=-gauss->weight*Jdettria*gamma*K[0]*dt; //don't forget the - + DL_scalar3=-gauss->weight*Jdettria*gamma*K[1]*dt; //don't forget the - + + } + else{ + DL_scalar2=-gauss->weight*Jdettria*gamma*K[0]; //don't forget the - + DL_scalar3=-gauss->weight*Jdettria*gamma*K[1]; //don't forget the - + } + + TripleMultiply( &dLx[0],1,numdof,1, &DL_scalar2,1,1,0, &L[0],1,numdof,0, &Ke_gg_gaussian2[0][0],0); + TripleMultiply( &dLy[0],1,numdof,1, &DL_scalar3,1,1,0, &L[0],1,numdof,0, &Ke_gg_gaussian3[0][0],0); + + for(i=0;ivalues[i*numdof+j]+=Ke_gg_gaussian1[i][j]; + for(i=0;ivalues[i*numdof+j]+=Ke_gg_gaussian2[i][j]; + for(i=0;ivalues[i*numdof+j]+=Ke_gg_gaussian3[i][j]; + } + + /*Clean up and return*/ + delete gauss; + return Ke; +} +/*}}}*/ /*FUNCTION Tria::CreateKMatrixPrognostic {{{1*/ ElementMatrix* Tria::CreateKMatrixPrognostic(void){ @@ -1530,4 +1614,7 @@ case PrognosticAnalysisEnum: pe=CreatePVectorPrognostic(); + break; + case HydrologyAnalysisEnum: + pe=CreatePVectorHydrology(); break; case BalancedthicknessAnalysisEnum: @@ -2293,4 +2380,53 @@ } /*}}}*/ +/*FUNCTION Tria::CreatePVectorHydrology {{{1*/ +ElementVector* Tria::CreatePVectorHydrology(void){ + + /*Constants*/ + const int numdof=NDOF1*NUMVERTICES; + + /*Intermediaries */ + int i,j,ig; + double Jdettria,dt; + double melting_g; + double old_watercolumn_g; + double xyz_list[NUMVERTICES][3]; + double L[NUMVERTICES]; + GaussTria* gauss=NULL; + + /*Initialize Element vector*/ + ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + this->parameters->FindParam(&dt,DtEnum); + Input* melting_input=inputs->GetInput(MeltingRateEnum); _assert_(melting_input); + Input* old_watercolumn_input=inputs->GetInput(OldWaterColumnEnum); _assert_(old_watercolumn_input); + + /*Initialize melting_correction_g to 0, do not forget!:*/ + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(2); + for(ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); + GetL(&L[0], &xyz_list[0][0], gauss,NDOF1); + + melting_input->GetParameterValue(&melting_g,gauss); + old_watercolumn_input->GetParameterValue(&old_watercolumn_g,gauss); + + + + if(dt)for(i=0;ivalues[i]+=Jdettria*gauss->weight*(old_watercolumn_g+dt*melting_g)*L[i]; + else for(i=0;ivalues[i]+=Jdettria*gauss->weight*melting_g*L[i]; + + } + + /*Clean up and return*/ + delete gauss; + return pe; +} +/*}}}*/ /*FUNCTION Tria::CreatePVectorPrognostic_CG {{{1*/ ElementVector* Tria::CreatePVectorPrognostic_CG(void){ @@ -2895,4 +3031,6 @@ else if (analysis_type==DiagnosticHutterAnalysisEnum) GetSolutionFromInputsDiagnosticHutter(solution); + else if (analysis_type==HydrologyAnalysisEnum) + GetSolutionFromInputsHydrology(solution); else _error_("analysis: %s not supported yet",EnumToString(analysis_type)); @@ -2969,4 +3107,41 @@ values[i*NDOF2+0]=vx; values[i*NDOF2+1]=vy; + } + + VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); + + /*Free ressources:*/ + delete gauss; + xfree((void**)&doflist); +} +/*}}}*/ +/*FUNCTION Tria::GetSolutionFromInputsHydrology{{{1*/ +void Tria::GetSolutionFromInputsHydrology(Vec solution){ + + const int numdof=NDOF1*NUMVERTICES; + + int i,dummy; + int* doflist=NULL; + double watercolumn; + double values[numdof]; + GaussTria* gauss=NULL; + + /*Get dof list: */ + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + + /*Get inputs*/ + Input* watercolumn_input=inputs->GetInput(WaterColumnEnum); _assert_(watercolumn_input); + + /*Ok, we have watercolumn values, fill in watercolumn array: */ + /*P1 element only for now*/ + gauss=new GaussTria(); + for(i=0;iGaussVertex(i); + + /*Recover watercolumn*/ + watercolumn_input->GetParameterValue(&watercolumn,gauss); + values[i]=watercolumn; + if((values[i])<0)values[i]=0; } @@ -3015,4 +3190,54 @@ } } +} +/*}}}*/ +/*FUNCTION Tria::GetHydrologyK {{{1*/ +void Tria::GetHydrologyK(double* K,double* xyz_list,GaussTria* gauss){ + + + /*material parameters: */ + double rho_ice; + double rho_water; + double g; + double dsdx, dsdy; + double dbdx, dbdy; + double surface_slope; + double kn; + double w; + + Input* surfaceslopex_input=NULL; + Input* surfaceslopey_input=NULL; + Input* bedslopex_input=NULL; + Input* bedslopey_input=NULL; + Input* watercolumn_input=NULL; + + + /*recover parameters and inputs: */ + rho_ice=matpar->GetRhoIce(); + rho_water=matpar->GetRhoWater(); + g=matpar->GetG(); + kn=matpar->GetKn(); + + surfaceslopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input); + surfaceslopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input); + bedslopex_input=inputs->GetInput(BedSlopeXEnum); _assert_(bedslopex_input); + bedslopey_input=inputs->GetInput(BedSlopeYEnum); _assert_(bedslopey_input); + watercolumn_input=inputs->GetInput(WaterColumnEnum); _assert_(watercolumn_input); + + /*recover slopes: */ + surfaceslopex_input->GetParameterValue(&dsdx,gauss); + surfaceslopey_input->GetParameterValue(&dsdy,gauss); + bedslopex_input->GetParameterValue(&dbdx,gauss); + bedslopey_input->GetParameterValue(&dbdy,gauss); + + /*magnitude of surface slope: */ + surface_slope=sqrt( pow(dsdx,2) + pow(dsdy,2)); + + /*recover water column: */ + watercolumn_input->GetParameterValue(&w,gauss); + + K[0]=pow(w,2)*(rho_ice*g*dsdx+(rho_water/rho_ice-1)*rho_ice*g*dbdx) - rho_ice * g * kn* w * (dsdx - dbdx ) * surface_slope; + K[1]=pow(w,2)*(rho_ice*g*dsdy+(rho_water/rho_ice-1)*rho_ice*g*dbdy) - rho_ice * g * kn *w * (dsdy - dbdy ) * surface_slope; + } /*}}}*/ @@ -3662,4 +3887,10 @@ this->inputs->AddInput(new TriaVertexInput(MeltingRateCorrectionEnum,nodeinputs)); } + if (iomodel->watercolumn){ + for(i=0;i<3;i++)nodeinputs[i]=iomodel->watercolumn[tria_vertex_ids[i]-1]; + this->inputs->AddInput(new TriaVertexInput(WaterColumnEnum,nodeinputs)); + this->inputs->AddInput(new TriaVertexInput(OldWaterColumnEnum,nodeinputs)); + } + if (iomodel->accumulation_rate) { for(i=0;i<3;i++)nodeinputs[i]=iomodel->accumulation_rate[tria_vertex_ids[i]-1]/iomodel->yts; @@ -3792,4 +4023,7 @@ InputUpdateFromSolutionPrognostic(solution); break; + case HydrologyAnalysisEnum: + InputUpdateFromSolutionHydrology(solution); + break; case BalancedthicknessAnalysisEnum: InputUpdateFromSolutionOneDof(solution,ThicknessEnum); @@ -4100,4 +4334,35 @@ } /*}}}*/ +/*FUNCTION Tria::InputUpdateFromSolutionHydrology{{{1*/ +void Tria::InputUpdateFromSolutionHydrology(double* solution){ + + /*Intermediaries*/ + const int numdof = NDOF1*NUMVERTICES; + + int i; + int* doflist=NULL; + double values[numdof]; + + /*Get dof list: */ + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + + /*Use the dof list to index into the solution vector: */ + for(i=0;iinputs->ChangeEnum(WaterColumnEnum,PicardWaterColumnEnum); + + /*Add input to the element: */ + this->inputs->AddInput(new TriaVertexInput(WaterColumnEnum,values)); + + /*Free ressources:*/ + xfree((void**)&doflist); +} +/*}}}*/ /*FUNCTION Tria::InputUpdateFromVector(double* vector, int name, int type);{{{1*/ void Tria::InputUpdateFromVector(double* vector, int name, int type){ @@ -4233,4 +4498,5 @@ name==SurfaceSlopeYEnum || name==MeltingRateEnum || + name==WaterColumnEnum || name==AccumulationRateEnum || name==SurfaceAreaEnum|| Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 7639) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 7640) @@ -154,4 +154,5 @@ ElementMatrix* CreateKMatrixDiagnosticVertSurface(void); ElementMatrix* CreateKMatrixMelting(void); + ElementMatrix* CreateKMatrixHydrology(void); ElementMatrix* CreateKMatrixPrognostic(void); ElementMatrix* CreateKMatrixPrognostic_CG(void); @@ -169,4 +170,5 @@ ElementVector* CreatePVectorAdjointBalancedthickness(void); ElementVector* CreatePVectorDiagnosticHutter(void); + ElementVector* CreatePVectorHydrology(void); ElementVector* CreatePVectorPrognostic(void); ElementVector* CreatePVectorPrognostic_CG(void); @@ -185,4 +187,5 @@ void GetSolutionFromInputsDiagnosticHoriz(Vec solution); void GetSolutionFromInputsDiagnosticHutter(Vec solution); + void GetSolutionFromInputsHydrology(Vec solution); void GetStrainRate2d(double* epsilon,double* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input); void GradjDragStokes(Vec gradient); @@ -193,7 +196,9 @@ void InputUpdateFromSolutionOneDof(double* solution,int enum_type); void InputUpdateFromSolutionPrognostic(double* solution); + void InputUpdateFromSolutionHydrology(double* solution); bool IsInput(int name); void SetClone(int* minranks); void SurfaceNormal(double* surface_normal, double xyz_list[3][3]); + void GetHydrologyK(double* K,double* xyz_list,GaussTria* gauss); /*}}}*/ Index: /issm/trunk/src/c/objects/IoModel.cpp =================================================================== --- /issm/trunk/src/c/objects/IoModel.cpp (revision 7639) +++ /issm/trunk/src/c/objects/IoModel.cpp (revision 7640) @@ -80,8 +80,10 @@ xfree((void**)&this->spcvelocity); xfree((void**)&this->spcthickness); + xfree((void**)&this->spcwatercolumn); xfree((void**)&this->spctemperature); xfree((void**)&this->edges); xfree((void**)&this->geothermalflux); xfree((void**)&this->melting_rate); + xfree((void**)&this->watercolumn); xfree((void**)&this->melting_rate_correction); xfree((void**)&this->accumulation_rate); @@ -199,4 +201,6 @@ /*!Get thermal parameters: */ IoModelFetchData(&this->beta,iomodel_handle,"beta"); + IoModelFetchData(&this->gamma,iomodel_handle,"gamma"); + IoModelFetchData(&this->kn,iomodel_handle,"kn"); IoModelFetchData(&this->meltingpoint,iomodel_handle,"meltingpoint"); IoModelFetchData(&this->latentheat,iomodel_handle,"latentheat"); @@ -274,4 +278,5 @@ this->gl_melting_rate=0; this->melting_rate=NULL; + this->watercolumn=NULL; this->melting_rate_correction=NULL; this->melting_rate_correction_apply=0; @@ -306,4 +311,5 @@ this-> spctemperature=NULL; this-> spcthickness=NULL; + this-> spcwatercolumn=NULL; this->numberofedges=0; this->edges=NULL; Index: /issm/trunk/src/c/objects/IoModel.h =================================================================== --- /issm/trunk/src/c/objects/IoModel.h (revision 7639) +++ /issm/trunk/src/c/objects/IoModel.h (revision 7640) @@ -103,4 +103,5 @@ double* spctemperature; double* spcthickness; + double* spcwatercolumn; double* geothermalflux; int numberofedges; @@ -163,4 +164,6 @@ /*thermal parameters: */ double beta; + double gamma; + double kn; double meltingpoint; double latentheat; @@ -182,4 +185,5 @@ /*basal: */ double* melting_rate; + double* watercolumn; double gl_melting_rate; double* melting_rate_correction; Index: /issm/trunk/src/c/objects/Materials/Matpar.cpp =================================================================== --- /issm/trunk/src/c/objects/Materials/Matpar.cpp (revision 7639) +++ /issm/trunk/src/c/objects/Materials/Matpar.cpp (revision 7640) @@ -35,4 +35,6 @@ double matpar_thermal_exchange_velocity; double matpar_g; + double matpar_gamma; + double matpar_kn; matpar_g=iomodel->g; @@ -46,4 +48,6 @@ matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity; matpar_thermal_exchange_velocity=iomodel->thermal_exchange_velocity; + matpar_gamma=iomodel->gamma; + matpar_kn=iomodel->kn; this->mid=matpar_mid; @@ -58,4 +62,6 @@ this->thermal_exchange_velocity=matpar_thermal_exchange_velocity; this->g=matpar_g; + this->gamma=matpar_gamma; + this->kn=matpar_kn; } @@ -83,4 +89,6 @@ printf(" thermal_exchange_velocity: %g\n",thermal_exchange_velocity); printf(" g: %g\n",g); + printf(" gamma: %g\n",gamma); + printf(" kn: %g\n",kn); return; } @@ -101,4 +109,6 @@ printf(" thermal_exchange_velocity: %g\n",thermal_exchange_velocity); printf(" g: %g\n",g); + printf(" gamma: %g\n",gamma); + printf(" kn: %g\n",kn); return; } @@ -140,4 +150,6 @@ memcpy(marshalled_dataset,&thermal_exchange_velocity,sizeof(thermal_exchange_velocity));marshalled_dataset+=sizeof(thermal_exchange_velocity); memcpy(marshalled_dataset,&g,sizeof(g));marshalled_dataset+=sizeof(g); + memcpy(marshalled_dataset,&gamma,sizeof(gamma));marshalled_dataset+=sizeof(gamma); + memcpy(marshalled_dataset,&kn,sizeof(kn));marshalled_dataset+=sizeof(kn); *pmarshalled_dataset=marshalled_dataset; @@ -159,4 +171,6 @@ sizeof(thermal_exchange_velocity)+ sizeof(g)+ + sizeof(gamma)+ + sizeof(kn)+ sizeof(int); //sizeof(int) for enum type } @@ -184,4 +198,6 @@ memcpy(&thermal_exchange_velocity,marshalled_dataset,sizeof(thermal_exchange_velocity));marshalled_dataset+=sizeof(thermal_exchange_velocity); memcpy(&g,marshalled_dataset,sizeof(g));marshalled_dataset+=sizeof(g); + memcpy(&gamma,marshalled_dataset,sizeof(gamma));marshalled_dataset+=sizeof(gamma); + memcpy(&kn,marshalled_dataset,sizeof(kn));marshalled_dataset+=sizeof(kn); /*return: */ @@ -359,4 +375,12 @@ } /*}}}1*/ - - +/*FUNCTION Matpar::GetGamma {{{1*/ +double Matpar::GetGamma(){ + return gamma; +} +/*}}}1*/ +/*FUNCTION Matpar::GetKn {{{1*/ +double Matpar::GetKn(){ + return kn; +} +/*}}}1*/ Index: /issm/trunk/src/c/objects/Materials/Matpar.h =================================================================== --- /issm/trunk/src/c/objects/Materials/Matpar.h (revision 7639) +++ /issm/trunk/src/c/objects/Materials/Matpar.h (revision 7640) @@ -26,4 +26,7 @@ double thermal_exchange_velocity; double g; + /*hydrology: */ + double kn; + double gamma; public: @@ -73,4 +76,6 @@ double GetBeta(); double GetMeltingPoint(); + double GetGamma(); + double GetKn(); /*}}}*/ Index: /issm/trunk/src/c/solutions/CorePointerFromSolutionEnum.cpp =================================================================== --- /issm/trunk/src/c/solutions/CorePointerFromSolutionEnum.cpp (revision 7639) +++ /issm/trunk/src/c/solutions/CorePointerFromSolutionEnum.cpp (revision 7640) @@ -58,5 +58,7 @@ solutioncore=&groundinglinemigration2d_core; break; - + case HydrologySolutionEnum: + solutioncore=&hydrology_core; + break; default: _error_("%s%s%s"," solution type: ",EnumToString(solutiontype)," not supported yet!"); Index: /issm/trunk/src/c/solutions/SolutionConfiguration.cpp =================================================================== --- /issm/trunk/src/c/solutions/SolutionConfiguration.cpp (revision 7639) +++ /issm/trunk/src/c/solutions/SolutionConfiguration.cpp (revision 7640) @@ -59,4 +59,12 @@ analyses[0]=ThermalAnalysisEnum; analyses[1]=MeltingAnalysisEnum; + break; + + case HydrologySolutionEnum: + numanalyses=3; + analyses=(int*)xmalloc(numanalyses*sizeof(int)); + analyses[0]=HydrologyAnalysisEnum; + analyses[1]=SurfaceSlopeAnalysisEnum; + analyses[2]=BedSlopeAnalysisEnum; break; Index: /issm/trunk/src/c/solutions/solutions.h =================================================================== --- /issm/trunk/src/c/solutions/solutions.h (revision 7639) +++ /issm/trunk/src/c/solutions/solutions.h (revision 7640) @@ -17,4 +17,6 @@ void gradient_core(FemModel* femmodel,int step=0, double search_scalar=0); void diagnostic_core(FemModel* femmodel); +void hydrology_core(FemModel* femmodel); +void hydrology_core_step(FemModel* femmodel,int step, double time); void thermal_core(FemModel* femmodel); void thermal_core_step(FemModel* femmodel,int step, double time);