Index: /issm/trunk-jpl/src/c/classes/Constraints/SpcDynamic.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Constraints/SpcDynamic.cpp (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Constraints/SpcDynamic.cpp (revision 20810) @@ -71,10 +71,4 @@ int SpcDynamic::Id(void){ return sid; }/*{{{*/ /*}}}*/ -int SpcDynamic::ObjectEnum(void){/*{{{*/ - - return SpcDynamicEnum; - -} -/*}}}*/ void SpcDynamic::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/ @@ -88,4 +82,10 @@ MARSHALLING(isset); MARSHALLING(penalty); + +} +/*}}}*/ +int SpcDynamic::ObjectEnum(void){/*{{{*/ + + return SpcDynamicEnum; } Index: /issm/trunk-jpl/src/c/classes/Constraints/SpcDynamic.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Constraints/SpcDynamic.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Constraints/SpcDynamic.h (revision 20810) @@ -34,6 +34,6 @@ void Echo(); int Id(); + void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction); int ObjectEnum(); - void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction); /*Constraint virtual functions definitions*/ Index: /issm/trunk-jpl/src/c/classes/Constraints/SpcStatic.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Constraints/SpcStatic.cpp (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Constraints/SpcStatic.cpp (revision 20810) @@ -73,10 +73,4 @@ int SpcStatic::Id(void){ return sid; }/*{{{*/ /*}}}*/ -int SpcStatic::ObjectEnum(void){/*{{{*/ - - return SpcStaticEnum; - -} -/*}}}*/ void SpcStatic::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/ @@ -89,4 +83,10 @@ MARSHALLING(analysis_type); MARSHALLING(penalty); + +} +/*}}}*/ +int SpcStatic::ObjectEnum(void){/*{{{*/ + + return SpcStaticEnum; } Index: /issm/trunk-jpl/src/c/classes/Constraints/SpcStatic.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Constraints/SpcStatic.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Constraints/SpcStatic.h (revision 20810) @@ -29,10 +29,10 @@ /*}}}*/ /*Object virtual functions definitions:{{{ */ + Object* copy(); + void DeepEcho(); void Echo(); - void DeepEcho(); int Id(); + void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction); int ObjectEnum(); - Object* copy(); - void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction); /*}}}*/ /*Constraint virtual functions definitions: {{{*/ @@ -43,6 +43,6 @@ /*}}}*/ /*SpcStatic management:{{{ */ + int GetDof(); int GetNodeId(); - int GetDof(); IssmDouble GetValue(); /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Constraints/SpcTransient.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Constraints/SpcTransient.cpp (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Constraints/SpcTransient.cpp (revision 20810) @@ -79,10 +79,4 @@ } /*}}}*/ -int SpcTransient::ObjectEnum(void){/*{{{*/ - - return SpcTransientEnum; - -} -/*}}}*/ void SpcTransient::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/ @@ -106,4 +100,10 @@ } /*}}}*/ +int SpcTransient::ObjectEnum(void){/*{{{*/ + + return SpcTransientEnum; + +} +/*}}}*/ /*Constraint virtual functions definitions:*/ Index: /issm/trunk-jpl/src/c/classes/Constraints/SpcTransient.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Constraints/SpcTransient.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Constraints/SpcTransient.h (revision 20810) @@ -35,6 +35,6 @@ void Echo(); int Id(); + void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction); int ObjectEnum(); - void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction); /*}}}*/ /*Constraint virtual functions definitions: {{{*/ Index: /issm/trunk-jpl/src/c/classes/Contour.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Contour.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Contour.h (revision 20810) @@ -52,4 +52,15 @@ /*Object virtual function resolutoin: */ + Object* copy() {/*{{{*/ + + Contour* contour = new Contour(this->id,this->nods,this->x,this->y,this->closed); + + return (Object*) contour; + } + /*}}}*/ + void DeepEcho(void){/*{{{*/ + this->Echo(); + } + /*}}}*/ void Echo(void){/*{{{*/ _printf_(" Contour: " << id << "\n"); @@ -64,21 +75,6 @@ } /*}}}*/ - void DeepEcho(void){/*{{{*/ - this->Echo(); - } - /*}}}*/ int Id(void){/*{{{*/ return id; - } - /*}}}*/ - int ObjectEnum(void){/*{{{*/ - return ContourEnum; - } - /*}}}*/ - Object* copy() {/*{{{*/ - - Contour* contour = new Contour(this->id,this->nods,this->x,this->y,this->closed); - - return (Object*) contour; } /*}}}*/ @@ -87,4 +83,8 @@ } /*}}}*/ + int ObjectEnum(void){/*{{{*/ + return ContourEnum; + } + /*}}}*/ }; Index: /issm/trunk-jpl/src/c/classes/DependentObject.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/DependentObject.cpp (revision 20809) +++ /issm/trunk-jpl/src/c/classes/DependentObject.cpp (revision 20810) @@ -36,4 +36,11 @@ /*Object virtual functions definitions:*/ +Object* DependentObject::copy(void) { /*{{{*/ + return new DependentObject(name,type,index); +} /*}}}*/ +void DependentObject::DeepEcho(void){/*{{{*/ + this->Echo(); +} +/*}}}*/ void DependentObject::Echo(void){/*{{{*/ @@ -49,8 +56,4 @@ } /*}}}*/ -void DependentObject::DeepEcho(void){/*{{{*/ - this->Echo(); -} -/*}}}*/ int DependentObject::Id(void){ return -1; }/*{{{*/ /*}}}*/ @@ -61,7 +64,4 @@ } /*}}}*/ -Object* DependentObject::copy(void) { /*{{{*/ - return new DependentObject(name,type,index); -} /*}}}*/ /*DependentObject methods: */ Index: /issm/trunk-jpl/src/c/classes/DependentObject.h =================================================================== --- /issm/trunk-jpl/src/c/classes/DependentObject.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/DependentObject.h (revision 20810) @@ -27,9 +27,9 @@ /*}}}*/ /*Object virtual functions definitions:{{{ */ + Object* copy(void); + void DeepEcho(); void Echo(); - void DeepEcho(); int Id(); int ObjectEnum(); - Object* copy(void); void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){/*{{{*/ _error_("not implemented yet!"); Index: /issm/trunk-jpl/src/c/classes/DofIndexing.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/DofIndexing.cpp (revision 20809) +++ /issm/trunk-jpl/src/c/classes/DofIndexing.cpp (revision 20810) @@ -214,5 +214,32 @@ /*Some of the Object functionality: */ -void DofIndexing::Echo(void){/*{{{*/ +void DofIndexing::Activate(void){/*{{{*/ + + this->active = true; + + /*Constrain to 0. at this point*/ + for(int i=0;igsize;i++){ + this->f_set[i] = true; + this->s_set[i] = false; + this->svalues[i] = 0.; + } + return; +} +/*}}}*/ +void DofIndexing::Deactivate(void){/*{{{*/ + this->active = false; + + /*Constrain to 0. at this point*/ + for(int i=0;igsize;i++){ + this->f_set[i] = false; + this->s_set[i] = true; + this->svalues[i] = 0.; + } + return; +} +/*}}}*/ +void DofIndexing::DeepEcho(void){/*{{{*/ + + int i; _printf_("DofIndexing:\n"); @@ -223,9 +250,45 @@ _printf_(" active: " << active << "\n"); _printf_(" freeze: " << freeze << "\n"); -} -/*}}}*/ -void DofIndexing::DeepEcho(void){/*{{{*/ - - int i; + + _printf_(" set membership: f,s sets \n"); + for(i=0;issize << "): |"); + for(i=0;igsize;i++){ + if(this->s_set[i])_printf_(" " << svalues[i] << " |"); + } + _printf_("\n"); + + if(doftype){ + _printf_(" doftype: |"); + for(i=0;igsize << "): |"); + for(i=0;igsize;i++){ + _printf_(" " << gdoflist[i] << " |"); + } + _printf_("\n"); + + _printf_(" f_doflist (" << this->fsize << "): |"); + for(i=0;ifsize;i++){ + _printf_(" " << fdoflist[i] << " |"); + } + _printf_("\n"); + + _printf_(" s_doflist (" << this->ssize << "): |"); + for(i=0;issize;i++){ + _printf_(" " << sdoflist[i] << " |"); + } + _printf_("\n"); +} +/*}}}*/ +void DofIndexing::Echo(void){/*{{{*/ _printf_("DofIndexing:\n"); @@ -236,67 +299,4 @@ _printf_(" active: " << active << "\n"); _printf_(" freeze: " << freeze << "\n"); - - _printf_(" set membership: f,s sets \n"); - for(i=0;issize << "): |"); - for(i=0;igsize;i++){ - if(this->s_set[i])_printf_(" " << svalues[i] << " |"); - } - _printf_("\n"); - - if(doftype){ - _printf_(" doftype: |"); - for(i=0;igsize << "): |"); - for(i=0;igsize;i++){ - _printf_(" " << gdoflist[i] << " |"); - } - _printf_("\n"); - - _printf_(" f_doflist (" << this->fsize << "): |"); - for(i=0;ifsize;i++){ - _printf_(" " << fdoflist[i] << " |"); - } - _printf_("\n"); - - _printf_(" s_doflist (" << this->ssize << "): |"); - for(i=0;issize;i++){ - _printf_(" " << sdoflist[i] << " |"); - } - _printf_("\n"); -} -/*}}}*/ -void DofIndexing::Deactivate(void){/*{{{*/ - this->active = false; - - /*Constrain to 0. at this point*/ - for(int i=0;igsize;i++){ - this->f_set[i] = false; - this->s_set[i] = true; - this->svalues[i] = 0.; - } - return; -} -/*}}}*/ -void DofIndexing::Activate(void){/*{{{*/ - - this->active = true; - - /*Constrain to 0. at this point*/ - for(int i=0;igsize;i++){ - this->f_set[i] = true; - this->s_set[i] = false; - this->svalues[i] = 0.; - } - return; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/DofIndexing.h =================================================================== --- /issm/trunk-jpl/src/c/classes/DofIndexing.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/DofIndexing.h (revision 20810) @@ -45,13 +45,13 @@ /*}}}*/ /*Object like functionality: {{{*/ + void copy(const DofIndexing& in); + void DeepEcho(void); void Echo(void); - void DeepEcho(void); - void copy(const DofIndexing& in); void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction); /*}}}*/ /*DofIndexing management: {{{*/ - DofIndexing* Spawn(int* indices, int numindices); void Activate(void); void Deactivate(void); + DofIndexing* Spawn(int* indices, int numindices); /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 20810) @@ -43,4 +43,92 @@ this->inputs->AddInput(input_in); }/*}}}*/ +void Element::ComputeLambdaS(){/*{{{*/ + + /*Intermediaries*/ + IssmDouble vx,vy,vz,vmag; + IssmDouble dvx[3],dvy[3],dvz[3],dvmag[3]; + IssmDouble epso,epsprime; + int dim; + IssmDouble *xyz_list = NULL; + + /*Retrieve all inputs we will be needing: */ + this->GetVerticesCoordinates(&xyz_list); + parameters->FindParam(&dim,DomainDimensionEnum); + Input* vx_input=this->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=this->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input=NULL; + if(dim==3){vz_input=this->GetInput(VzEnum); _assert_(vz_input);} + + /*Allocate arrays*/ + int numvertices = this->GetNumberOfVertices(); + IssmDouble* lambdas = xNew(numvertices); + IssmDouble* vorticityx = xNew(numvertices); + IssmDouble* vorticityy = xNew(numvertices); + IssmDouble* vorticityz = xNew(numvertices); + IssmDouble* omega_corrx = xNew(numvertices); + IssmDouble* omega_corry = xNew(numvertices); + IssmDouble* omega_corrz = xNew(numvertices); + + /* Start looping on the number of vertices: */ + Gauss* gauss=this->NewGauss(); + for (int iv=0;ivGaussVertex(iv); + + /*Get velocity derivatives in all directions*/ + _assert_(dim>1); + _assert_(vx_input); + vx_input->GetInputValue(&vx,gauss); + vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + _assert_(vy_input); + vy_input->GetInputValue(&vy,gauss); + vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + if(dim==3){ + _assert_(vz_input); + vz_input->GetInputValue(&vz,gauss); + vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss); + } + else{ + vz = 0.; + dvz[0] = 0.; dvz[1] = 0.; dvz[2] = 0.; + } + /*Calculate velocity magnitude and its derivative*/ + vmag = sqrt(vx*vx+vy*vy+vz*vz); + if(vmag<1e-12){ + vmag=1e-12; + dvmag[0]=0; + dvmag[1]=0; + dvmag[2]=0; + } + else{ + dvmag[0]=1./(2*sqrt(vmag))*(2*vx*dvx[0]+2*vy*dvy[0]+2*vz*dvz[0]); + dvmag[1]=1./(2*sqrt(vmag))*(2*vx*dvx[1]+2*vy*dvy[1]+2*vz*dvz[1]); + dvmag[2]=1./(2*sqrt(vmag))*(2*vx*dvx[2]+2*vy*dvy[2]+2*vz*dvz[2]); + } + EstarStrainrateQuantities(&epso,&epsprime,vx,vy,vz,vmag,&dvx[0],&dvy[0],&dvz[0],&dvmag[0]); + lambdas[iv]=EstarLambdaS(epso,epsprime); + vorticityx[iv]=dvz[1]-dvy[2]; + vorticityy[iv]=dvx[2]-dvz[0]; + vorticityz[iv]=dvy[0]-dvx[1]; + omega_corrx[iv] = -vorticityx[iv] + 2*(dvmag[1]*vz - dvmag[2]*vy)/vmag; + omega_corry[iv] = -vorticityy[iv] + 2*(dvmag[2]*vx - dvmag[0]*vz)/vmag; + omega_corrz[iv] = -vorticityz[iv] + 2*(dvmag[0]*vy - dvmag[1]*vx)/vmag; + } + + /*Add Stress tensor components into inputs*/ + this->AddInput(LambdaSEnum,lambdas,P1Enum); + this->AddInput(Outputdefinition21Enum,vorticityx,P1Enum); + this->AddInput(Outputdefinition22Enum,vorticityy,P1Enum); + this->AddInput(Outputdefinition23Enum,vorticityz,P1Enum); + this->AddInput(Outputdefinition31Enum,omega_corrx,P1Enum); + this->AddInput(Outputdefinition32Enum,omega_corry,P1Enum); + this->AddInput(Outputdefinition33Enum,omega_corrz,P1Enum); + + /*Clean up and return*/ + delete gauss; + xDelete(xyz_list); + xDelete(lambdas); + +} +/*}}}*/ void Element::ComputeNewDamage(){/*{{{*/ @@ -218,92 +306,4 @@ xDelete(eps_yz); xDelete(eps_ef); - -} -/*}}}*/ -void Element::ComputeLambdaS(){/*{{{*/ - - /*Intermediaries*/ - IssmDouble vx,vy,vz,vmag; - IssmDouble dvx[3],dvy[3],dvz[3],dvmag[3]; - IssmDouble epso,epsprime; - int dim; - IssmDouble *xyz_list = NULL; - - /*Retrieve all inputs we will be needing: */ - this->GetVerticesCoordinates(&xyz_list); - parameters->FindParam(&dim,DomainDimensionEnum); - Input* vx_input=this->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=this->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input=NULL; - if(dim==3){vz_input=this->GetInput(VzEnum); _assert_(vz_input);} - - /*Allocate arrays*/ - int numvertices = this->GetNumberOfVertices(); - IssmDouble* lambdas = xNew(numvertices); - IssmDouble* vorticityx = xNew(numvertices); - IssmDouble* vorticityy = xNew(numvertices); - IssmDouble* vorticityz = xNew(numvertices); - IssmDouble* omega_corrx = xNew(numvertices); - IssmDouble* omega_corry = xNew(numvertices); - IssmDouble* omega_corrz = xNew(numvertices); - - /* Start looping on the number of vertices: */ - Gauss* gauss=this->NewGauss(); - for (int iv=0;ivGaussVertex(iv); - - /*Get velocity derivatives in all directions*/ - _assert_(dim>1); - _assert_(vx_input); - vx_input->GetInputValue(&vx,gauss); - vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); - _assert_(vy_input); - vy_input->GetInputValue(&vy,gauss); - vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); - if(dim==3){ - _assert_(vz_input); - vz_input->GetInputValue(&vz,gauss); - vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss); - } - else{ - vz = 0.; - dvz[0] = 0.; dvz[1] = 0.; dvz[2] = 0.; - } - /*Calculate velocity magnitude and its derivative*/ - vmag = sqrt(vx*vx+vy*vy+vz*vz); - if(vmag<1e-12){ - vmag=1e-12; - dvmag[0]=0; - dvmag[1]=0; - dvmag[2]=0; - } - else{ - dvmag[0]=1./(2*sqrt(vmag))*(2*vx*dvx[0]+2*vy*dvy[0]+2*vz*dvz[0]); - dvmag[1]=1./(2*sqrt(vmag))*(2*vx*dvx[1]+2*vy*dvy[1]+2*vz*dvz[1]); - dvmag[2]=1./(2*sqrt(vmag))*(2*vx*dvx[2]+2*vy*dvy[2]+2*vz*dvz[2]); - } - EstarStrainrateQuantities(&epso,&epsprime,vx,vy,vz,vmag,&dvx[0],&dvy[0],&dvz[0],&dvmag[0]); - lambdas[iv]=EstarLambdaS(epso,epsprime); - vorticityx[iv]=dvz[1]-dvy[2]; - vorticityy[iv]=dvx[2]-dvz[0]; - vorticityz[iv]=dvy[0]-dvx[1]; - omega_corrx[iv] = -vorticityx[iv] + 2*(dvmag[1]*vz - dvmag[2]*vy)/vmag; - omega_corry[iv] = -vorticityy[iv] + 2*(dvmag[2]*vx - dvmag[0]*vz)/vmag; - omega_corrz[iv] = -vorticityz[iv] + 2*(dvmag[0]*vy - dvmag[1]*vx)/vmag; - } - - /*Add Stress tensor components into inputs*/ - this->AddInput(LambdaSEnum,lambdas,P1Enum); - this->AddInput(Outputdefinition21Enum,vorticityx,P1Enum); - this->AddInput(Outputdefinition22Enum,vorticityy,P1Enum); - this->AddInput(Outputdefinition23Enum,vorticityz,P1Enum); - this->AddInput(Outputdefinition31Enum,omega_corrx,P1Enum); - this->AddInput(Outputdefinition32Enum,omega_corry,P1Enum); - this->AddInput(Outputdefinition33Enum,omega_corrz,P1Enum); - - /*Clean up and return*/ - delete gauss; - xDelete(xyz_list); - xDelete(lambdas); } @@ -426,155 +426,4 @@ delete this->material; }/*}}}*/ -IssmDouble Element::Divergence(void){/*{{{*/ - /*Compute element divergence*/ - - /*Intermediaries*/ - int dim; - IssmDouble Jdet; - IssmDouble divergence=0.; - IssmDouble dvx[3],dvy[3],dvz[3]; - IssmDouble *xyz_list = NULL; - - /*Get inputs and parameters*/ - this->FindParam(&dim,DomainDimensionEnum); - Input* vx_input = this->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input = this->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input = NULL; - if(dim==3){ - vz_input = this->GetInput(VzEnum); _assert_(vz_input); - } - this->GetVerticesCoordinates(&xyz_list); - - Gauss* gauss=this->NewGauss(5); - for(int ig=gauss->begin();igend();ig++){ - gauss->GaussPoint(ig); - this->JacobianDeterminant(&Jdet,xyz_list,gauss); - - /*Get strain rate assuming that epsilon has been allocated*/ - vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); - vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); - if(dim==2){ - divergence += (dvx[0]+dvy[1])*gauss->weight*Jdet; - } - else{ - vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss); - divergence += (dvx[0]+dvy[1]+dvz[2])*gauss->weight*Jdet; - } - - } - - /*Clean up and return*/ - xDelete(xyz_list); - delete gauss; - return divergence; -}/*}}}*/ -void Element::dViscositydBFS(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/ - - /*Intermediaries*/ - IssmDouble dmudB; - IssmDouble epsilon3d[6];/* epsilon=[exx,eyy,exy,exy,exz,eyz]; */ - IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy]; */ - IssmDouble eps_eff; - IssmDouble eps0=1.e-27; - - if(dim==3){ - /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */ - this->StrainRateFS(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input,vz_input); - eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[3]*epsilon3d[3] + epsilon3d[4]*epsilon3d[4] + epsilon3d[5]*epsilon3d[5] + epsilon3d[0]*epsilon3d[1]+eps0*eps0); - } - else{ - /* eps_eff^2 = 1/2 ( exx^2 + eyy^2 + 2*exy^2 )*/ - this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); - eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + 2.*epsilon2d[2]*epsilon2d[2]); - } - /*Get viscosity*/ - material->GetViscosity_B(&dmudB,eps_eff); - - /*Assign output pointer*/ - *pdmudB=dmudB; - -} -/*}}}*/ -void Element::dViscositydBHO(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ - - /*Intermediaries*/ - IssmDouble dmudB; - IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exy,exz,eyz]; */ - IssmDouble epsilon2d[2];/* epsilon=[exx,eyy,exy]; */ - IssmDouble eps_eff; - IssmDouble eps0=1.e-27; - - if(dim==3){ - /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */ - this->StrainRateHO(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input); - eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[2]*epsilon3d[2] + epsilon3d[3]*epsilon3d[3] + epsilon3d[4]*epsilon3d[4] + epsilon3d[0]*epsilon3d[1]+eps0*eps0); - } - else{ - /* eps_eff^2 = 1/2 ( exx^2 + eyy^2 + 2*exy^2 )*/ - this->StrainRateHO2dvertical(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); - eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + 2.*epsilon2d[1]*epsilon2d[1] + eps0*eps0); - } - /*Get viscosity*/ - material->GetViscosity_B(&dmudB,eps_eff); - - /*Assign output pointer*/ - *pdmudB=dmudB; - -} -/*}}}*/ -void Element::dViscositydBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ - - /*Intermediaries*/ - IssmDouble dmudB; - IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy]; */ - IssmDouble epsilon1d; /* epsilon=[exx]; */ - IssmDouble eps_eff; - - if(dim==2){ - /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/ - this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); - eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]); - } - else{ - /* eps_eff^2 = 1/2 exx^2*/ - this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input); - eps_eff = sqrt(epsilon1d*epsilon1d/2.); - } - - /*Get viscosity*/ - material->GetViscosity_B(&dmudB,eps_eff); - - /*Assign output pointer*/ - *pdmudB=dmudB; - -} -/*}}}*/ -void Element::dViscositydDSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ - - /*Intermediaries*/ - IssmDouble dmudB; - IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy]; */ - IssmDouble epsilon1d; /* epsilon=[exx]; */ - IssmDouble eps_eff; - - if(dim==2){ - /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/ - this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); - eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]); - } - else{ - /* eps_eff^2 = 1/2 exx^2*/ - this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input); - eps_eff = sqrt(epsilon1d*epsilon1d/2.); - } - - /*Get viscosity*/ - material->GetViscosity_D(&dmudB,eps_eff); - - /*Assign output pointer*/ - *pdmudB=dmudB; - -} -/*}}}*/ void Element::Delta18oParameterization(void){/*{{{*/ @@ -681,103 +530,4 @@ } /*}}}*/ -void Element::MungsmtpParameterization(void){/*{{{*/ - /*Are we on the base? If not, return*/ - if(!IsOnBase()) return; - - int numvertices = this->GetNumberOfVertices(); - - int i; - IssmDouble* monthlytemperatures=xNew(12*numvertices); - IssmDouble* monthlyprec=xNew(12*numvertices); - IssmDouble* TemperaturesPresentday=xNew(12*numvertices); - IssmDouble* TemperaturesLgm=xNew(12*numvertices); - IssmDouble* PrecipitationsPresentday=xNew(12*numvertices); - IssmDouble* PrecipitationsLgm=xNew(12*numvertices); - IssmDouble* tmp=xNew(numvertices); - IssmDouble TdiffTime,PfacTime; - IssmDouble time,yts,time_yr; - this->parameters->FindParam(&time,TimeEnum); - this->parameters->FindParam(&yts,ConstantsYtsEnum); - time_yr=floor(time/yts)*yts; - - /*Recover present day temperature and precipitation*/ - Input* input=this->inputs->GetInput(SmbTemperaturesPresentdayEnum); _assert_(input); - Input* input2=this->inputs->GetInput(SmbTemperaturesLgmEnum); _assert_(input2); - Input* input3=this->inputs->GetInput(SmbPrecipitationsPresentdayEnum); _assert_(input3); - Input* input4=this->inputs->GetInput(SmbPrecipitationsLgmEnum); _assert_(input4); - /*loop over vertices: */ - Gauss* gauss=this->NewGauss(); - for(int month=0;month<12;month++) { - for(int iv=0;ivGaussVertex(iv); - input->GetInputValue(&TemperaturesPresentday[iv*12+month],gauss,month/12.*yts); - input2->GetInputValue(&TemperaturesLgm[iv*12+month],gauss,month/12.*yts); - input3->GetInputValue(&PrecipitationsPresentday[iv*12+month],gauss,month/12.*yts); - input4->GetInputValue(&PrecipitationsLgm[iv*12+month],gauss,month/12.*yts); - - PrecipitationsPresentday[iv*12+month]=PrecipitationsPresentday[iv*12+month]*yts; - PrecipitationsLgm[iv*12+month]=PrecipitationsLgm[iv*12+month]*yts; - } - } - - /*Recover interpolation parameters at time t*/ - this->parameters->FindParam(&TdiffTime,SmbTdiffEnum,time); - this->parameters->FindParam(&PfacTime,SmbPfacEnum,time); - - /*Compute the temperature and precipitation*/ - for(int iv=0;ivObjectEnum()){ - case TriaEnum: NewTemperatureInput->AddTimeInput(new TriaInput(SmbMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; - case PentaEnum: NewTemperatureInput->AddTimeInput(new PentaInput(SmbMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; - case TetraEnum: NewTemperatureInput->AddTimeInput(new TetraInput(SmbMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; - default: _error_("Not implemented yet"); - } - for(i=0;iObjectEnum()){ - case TriaEnum: NewPrecipitationInput->AddTimeInput(new TriaInput(SmbPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; - case PentaEnum: NewPrecipitationInput->AddTimeInput(new PentaInput(SmbPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; - case TetraEnum: NewPrecipitationInput->AddTimeInput(new TetraInput(SmbPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; - default: _error_("Not implemented yet"); - } - } - NewTemperatureInput->Configure(this->parameters); - NewPrecipitationInput->Configure(this->parameters); - - this->inputs->AddInput(NewTemperatureInput); - this->inputs->AddInput(NewPrecipitationInput); - - switch(this->ObjectEnum()){ - case TriaEnum: break; - case PentaEnum: - case TetraEnum: - this->InputExtrude(SmbMonthlytemperaturesEnum,-1); - this->InputExtrude(SmbPrecipitationEnum,-1); - break; - default: _error_("Not implemented yet"); - } - - /*clean-up*/ - delete gauss; - xDelete(monthlytemperatures); - xDelete(monthlyprec); - xDelete(TemperaturesPresentday); - xDelete(TemperaturesLgm); - xDelete(PrecipitationsPresentday); - xDelete(PrecipitationsLgm); - xDelete(tmp); - -} -/*}}}*/ void Element::Delta18opdParameterization(void){/*{{{*/ /*Are we on the base? If not, return*/ @@ -870,4 +620,155 @@ xDelete(tmp); +} +/*}}}*/ +IssmDouble Element::Divergence(void){/*{{{*/ + /*Compute element divergence*/ + + /*Intermediaries*/ + int dim; + IssmDouble Jdet; + IssmDouble divergence=0.; + IssmDouble dvx[3],dvy[3],dvz[3]; + IssmDouble *xyz_list = NULL; + + /*Get inputs and parameters*/ + this->FindParam(&dim,DomainDimensionEnum); + Input* vx_input = this->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = this->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input = NULL; + if(dim==3){ + vz_input = this->GetInput(VzEnum); _assert_(vz_input); + } + this->GetVerticesCoordinates(&xyz_list); + + Gauss* gauss=this->NewGauss(5); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + this->JacobianDeterminant(&Jdet,xyz_list,gauss); + + /*Get strain rate assuming that epsilon has been allocated*/ + vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + if(dim==2){ + divergence += (dvx[0]+dvy[1])*gauss->weight*Jdet; + } + else{ + vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss); + divergence += (dvx[0]+dvy[1]+dvz[2])*gauss->weight*Jdet; + } + + } + + /*Clean up and return*/ + xDelete(xyz_list); + delete gauss; + return divergence; +}/*}}}*/ +void Element::dViscositydBFS(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/ + + /*Intermediaries*/ + IssmDouble dmudB; + IssmDouble epsilon3d[6];/* epsilon=[exx,eyy,exy,exy,exz,eyz]; */ + IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy]; */ + IssmDouble eps_eff; + IssmDouble eps0=1.e-27; + + if(dim==3){ + /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */ + this->StrainRateFS(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input,vz_input); + eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[3]*epsilon3d[3] + epsilon3d[4]*epsilon3d[4] + epsilon3d[5]*epsilon3d[5] + epsilon3d[0]*epsilon3d[1]+eps0*eps0); + } + else{ + /* eps_eff^2 = 1/2 ( exx^2 + eyy^2 + 2*exy^2 )*/ + this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); + eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + 2.*epsilon2d[2]*epsilon2d[2]); + } + /*Get viscosity*/ + material->GetViscosity_B(&dmudB,eps_eff); + + /*Assign output pointer*/ + *pdmudB=dmudB; + +} +/*}}}*/ +void Element::dViscositydBHO(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ + + /*Intermediaries*/ + IssmDouble dmudB; + IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exy,exz,eyz]; */ + IssmDouble epsilon2d[2];/* epsilon=[exx,eyy,exy]; */ + IssmDouble eps_eff; + IssmDouble eps0=1.e-27; + + if(dim==3){ + /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */ + this->StrainRateHO(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input); + eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[2]*epsilon3d[2] + epsilon3d[3]*epsilon3d[3] + epsilon3d[4]*epsilon3d[4] + epsilon3d[0]*epsilon3d[1]+eps0*eps0); + } + else{ + /* eps_eff^2 = 1/2 ( exx^2 + eyy^2 + 2*exy^2 )*/ + this->StrainRateHO2dvertical(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); + eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + 2.*epsilon2d[1]*epsilon2d[1] + eps0*eps0); + } + /*Get viscosity*/ + material->GetViscosity_B(&dmudB,eps_eff); + + /*Assign output pointer*/ + *pdmudB=dmudB; + +} +/*}}}*/ +void Element::dViscositydBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ + + /*Intermediaries*/ + IssmDouble dmudB; + IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy]; */ + IssmDouble epsilon1d; /* epsilon=[exx]; */ + IssmDouble eps_eff; + + if(dim==2){ + /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/ + this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); + eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]); + } + else{ + /* eps_eff^2 = 1/2 exx^2*/ + this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input); + eps_eff = sqrt(epsilon1d*epsilon1d/2.); + } + + /*Get viscosity*/ + material->GetViscosity_B(&dmudB,eps_eff); + + /*Assign output pointer*/ + *pdmudB=dmudB; + +} +/*}}}*/ +void Element::dViscositydDSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ + + /*Intermediaries*/ + IssmDouble dmudB; + IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy]; */ + IssmDouble epsilon1d; /* epsilon=[exx]; */ + IssmDouble eps_eff; + + if(dim==2){ + /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/ + this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); + eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]); + } + else{ + /* eps_eff^2 = 1/2 exx^2*/ + this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input); + eps_eff = sqrt(epsilon1d*epsilon1d/2.); + } + + /*Get viscosity*/ + material->GetViscosity_D(&dmudB,eps_eff); + + /*Assign output pointer*/ + *pdmudB=dmudB; + } /*}}}*/ @@ -1554,12 +1455,4 @@ bool Element::IsIceInElement(){/*{{{*/ return (this->inputs->Min(MaskIceLevelsetEnum)<0.); -} -/*}}}*/ -bool Element::IsWaterInElement(){/*{{{*/ - return (this->inputs->Max(MaskOceanLevelsetEnum)>0.); -} -/*}}}*/ -bool Element::IsLandInElement(){/*{{{*/ - return (this->inputs->Max(MaskLandLevelsetEnum)>0.); } /*}}}*/ @@ -1650,4 +1543,12 @@ } /*}}}*/ +bool Element::IsLandInElement(){/*{{{*/ + return (this->inputs->Max(MaskLandLevelsetEnum)>0.); +} +/*}}}*/ +bool Element::IsWaterInElement(){/*{{{*/ + return (this->inputs->Max(MaskOceanLevelsetEnum)>0.); +} +/*}}}*/ void Element::LinearFloatingiceMeltingRate(){/*{{{*/ @@ -1860,4 +1761,103 @@ }/*}}}*/ +void Element::MungsmtpParameterization(void){/*{{{*/ + /*Are we on the base? If not, return*/ + if(!IsOnBase()) return; + + int numvertices = this->GetNumberOfVertices(); + + int i; + IssmDouble* monthlytemperatures=xNew(12*numvertices); + IssmDouble* monthlyprec=xNew(12*numvertices); + IssmDouble* TemperaturesPresentday=xNew(12*numvertices); + IssmDouble* TemperaturesLgm=xNew(12*numvertices); + IssmDouble* PrecipitationsPresentday=xNew(12*numvertices); + IssmDouble* PrecipitationsLgm=xNew(12*numvertices); + IssmDouble* tmp=xNew(numvertices); + IssmDouble TdiffTime,PfacTime; + IssmDouble time,yts,time_yr; + this->parameters->FindParam(&time,TimeEnum); + this->parameters->FindParam(&yts,ConstantsYtsEnum); + time_yr=floor(time/yts)*yts; + + /*Recover present day temperature and precipitation*/ + Input* input=this->inputs->GetInput(SmbTemperaturesPresentdayEnum); _assert_(input); + Input* input2=this->inputs->GetInput(SmbTemperaturesLgmEnum); _assert_(input2); + Input* input3=this->inputs->GetInput(SmbPrecipitationsPresentdayEnum); _assert_(input3); + Input* input4=this->inputs->GetInput(SmbPrecipitationsLgmEnum); _assert_(input4); + /*loop over vertices: */ + Gauss* gauss=this->NewGauss(); + for(int month=0;month<12;month++) { + for(int iv=0;ivGaussVertex(iv); + input->GetInputValue(&TemperaturesPresentday[iv*12+month],gauss,month/12.*yts); + input2->GetInputValue(&TemperaturesLgm[iv*12+month],gauss,month/12.*yts); + input3->GetInputValue(&PrecipitationsPresentday[iv*12+month],gauss,month/12.*yts); + input4->GetInputValue(&PrecipitationsLgm[iv*12+month],gauss,month/12.*yts); + + PrecipitationsPresentday[iv*12+month]=PrecipitationsPresentday[iv*12+month]*yts; + PrecipitationsLgm[iv*12+month]=PrecipitationsLgm[iv*12+month]*yts; + } + } + + /*Recover interpolation parameters at time t*/ + this->parameters->FindParam(&TdiffTime,SmbTdiffEnum,time); + this->parameters->FindParam(&PfacTime,SmbPfacEnum,time); + + /*Compute the temperature and precipitation*/ + for(int iv=0;ivObjectEnum()){ + case TriaEnum: NewTemperatureInput->AddTimeInput(new TriaInput(SmbMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; + case PentaEnum: NewTemperatureInput->AddTimeInput(new PentaInput(SmbMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; + case TetraEnum: NewTemperatureInput->AddTimeInput(new TetraInput(SmbMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; + default: _error_("Not implemented yet"); + } + for(i=0;iObjectEnum()){ + case TriaEnum: NewPrecipitationInput->AddTimeInput(new TriaInput(SmbPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; + case PentaEnum: NewPrecipitationInput->AddTimeInput(new PentaInput(SmbPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; + case TetraEnum: NewPrecipitationInput->AddTimeInput(new TetraInput(SmbPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; + default: _error_("Not implemented yet"); + } + } + NewTemperatureInput->Configure(this->parameters); + NewPrecipitationInput->Configure(this->parameters); + + this->inputs->AddInput(NewTemperatureInput); + this->inputs->AddInput(NewPrecipitationInput); + + switch(this->ObjectEnum()){ + case TriaEnum: break; + case PentaEnum: + case TetraEnum: + this->InputExtrude(SmbMonthlytemperaturesEnum,-1); + this->InputExtrude(SmbPrecipitationEnum,-1); + break; + default: _error_("Not implemented yet"); + } + + /*clean-up*/ + delete gauss; + xDelete(monthlytemperatures); + xDelete(monthlyprec); + xDelete(TemperaturesPresentday); + xDelete(TemperaturesLgm); + xDelete(PrecipitationsPresentday); + xDelete(PrecipitationsLgm); + xDelete(tmp); + +} +/*}}}*/ ElementMatrix* Element::NewElementMatrix(int approximation_enum){/*{{{*/ return new ElementMatrix(nodes,this->GetNumberOfNodes(),this->parameters,approximation_enum); @@ -2069,5 +2069,4 @@ switch(calvinglaw){ case DefaultCalvingEnum: - case CalvingMinthicknessEnum: //do nothing break; Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 20810) @@ -63,10 +63,8 @@ void ComputeStrainRate(); void CoordinateSystemTransform(IssmDouble** ptransform,Node** nodes,int numnodes,int* cs_array); - void Echo(); void DeepEcho(); void DeleteInput(int input_enum); void DeleteMaterials(void); void Delta18oParameterization(void); - void MungsmtpParameterization(void); void Delta18opdParameterization(void); IssmDouble Divergence(void); @@ -75,4 +73,5 @@ void dViscositydBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); void dViscositydDSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); + void Echo(); IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure); IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure); @@ -120,7 +119,7 @@ bool IsFloating(); bool IsIceInElement(); + bool IsInput(int name); + bool IsLandInElement(); bool IsWaterInElement(); - bool IsLandInElement(); - bool IsInput(int name); void LinearFloatingiceMeltingRate(); void MantlePlumeGeothermalFlux(); @@ -128,4 +127,5 @@ void MigrateGroundingLine(IssmDouble* sheet_ungrounding); void MismipFloatingiceMeltingRate(); + void MungsmtpParameterization(void); ElementMatrix* NewElementMatrix(int approximation_enum=NoneApproximationEnum); ElementMatrix* NewElementMatrixCoupling(int number_nodes,int approximation_enum=NoneApproximationEnum); @@ -174,8 +174,6 @@ virtual void AddInput(int input_enum, IssmDouble* values, int interpolation_enum)=0; virtual void AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part)=0; + virtual void CalvingRateDev(void){_error_("not implemented yet");}; virtual void CalvingRateLevermann(void)=0; - virtual void CalvingRateDev(void){_error_("not implemented yet");}; - virtual void WriteLevelsetSegment(DataSet* segments){_error_("not implemented yet");}; - virtual void ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments){_error_("not implemented yet");}; virtual IssmDouble CharacteristicLength(void)=0; virtual void ComputeBasalStress(Vector* sigma_b)=0; @@ -266,4 +264,5 @@ virtual void ResetFSBasalBoundaryCondition()=0; virtual void ResetHooks()=0; + virtual void ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments){_error_("not implemented yet");}; virtual void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index)=0; virtual void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters)=0; @@ -291,4 +290,5 @@ virtual void VerticalSegmentIndices(int** pindices,int* pnumseg)=0; virtual void ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0; + virtual void WriteLevelsetSegment(DataSet* segments){_error_("not implemented yet");}; virtual void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum)=0; @@ -297,9 +297,9 @@ #endif #ifdef _HAVE_SEALEVELRISE_ + virtual IssmDouble GetArea3D(void)=0; + virtual IssmDouble OceanAverage(IssmDouble* Sg)=0; + virtual IssmDouble OceanArea(void)=0; virtual void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=0; virtual void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=0; - virtual IssmDouble OceanAverage(IssmDouble* Sg)=0; - virtual IssmDouble OceanArea(void)=0; - virtual IssmDouble GetArea3D(void)=0; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Elements.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Elements.cpp (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Elements/Elements.cpp (revision 20810) @@ -48,30 +48,10 @@ } /*}}}*/ -void Elements::SetCurrentConfiguration(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters){/*{{{*/ +void Elements::InputDuplicate(int input_enum,int output_enum){/*{{{*/ - vector::iterator object; - Element* element=NULL; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - element=xDynamicCast((*object)); - element->SetCurrentConfiguration(elements,loads,nodes,materials,parameters); - + for(int i=0;iSize();i++){ + Element* element=xDynamicCast(this->GetObjectByOffset(i)); + element->InputDuplicate(input_enum,output_enum); } - -} -/*}}}*/ -void Elements::ResetHooks(){/*{{{*/ - - vector::iterator object; - Element* element=NULL; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - element=xDynamicCast((*object)); - element->ResetHooks(); - - } - } /*}}}*/ @@ -108,10 +88,30 @@ } /*}}}*/ -void Elements::InputDuplicate(int input_enum,int output_enum){/*{{{*/ +void Elements::ResetHooks(){/*{{{*/ - for(int i=0;iSize();i++){ - Element* element=xDynamicCast(this->GetObjectByOffset(i)); - element->InputDuplicate(input_enum,output_enum); + vector::iterator object; + Element* element=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + element=xDynamicCast((*object)); + element->ResetHooks(); + } + } /*}}}*/ +void Elements::SetCurrentConfiguration(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters){/*{{{*/ + + vector::iterator object; + Element* element=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + element=xDynamicCast((*object)); + element->SetCurrentConfiguration(elements,loads,nodes,materials,parameters); + + } + +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Elements.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Elements.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Elements/Elements.h (revision 20810) @@ -25,9 +25,9 @@ /*numerics*/ void Configure(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters); + void InputDuplicate(int input_enum,int output_enum); int MaxNumNodes(void); + int NumberOfElements(void); void SetCurrentConfiguration(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters); void ResetHooks(); - int NumberOfElements(void); - void InputDuplicate(int input_enum,int output_enum); }; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 20810) @@ -3386,4 +3386,42 @@ #ifdef _HAVE_DAKOTA_ +void Penta::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*{{{*/ + + int i,t,row; + IssmDouble time; + TransientInput *transientinput = NULL; + IssmDouble values[6]; + + /*Check that name is an element input*/ + if (!IsInput(name)) return; + + switch(type){ + + case VertexEnum: + /*Create transient input: */ + for(t=0;tvertices[i]->Sid(); + values[i]=matrix[ncols*row+t]; + } + + /*time:*/ + time=matrix[(nrows-1)*ncols+t]; + + if(t==0) transientinput=new TransientInput(name); + transientinput->AddTimeInput(new PentaInput(name,values,P1Enum),time); + transientinput->Configure(parameters); + } + this->inputs->AddInput(transientinput); + break; + + default: + _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet"); + } + +} +/*}}}*/ void Penta::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*{{{*/ @@ -3486,41 +3524,3 @@ } /*}}}*/ -void Penta::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*{{{*/ - - int i,t,row; - IssmDouble time; - TransientInput *transientinput = NULL; - IssmDouble values[6]; - - /*Check that name is an element input*/ - if (!IsInput(name)) return; - - switch(type){ - - case VertexEnum: - /*Create transient input: */ - for(t=0;tvertices[i]->Sid(); - values[i]=matrix[ncols*row+t]; - } - - /*time:*/ - time=matrix[(nrows-1)*ncols+t]; - - if(t==0) transientinput=new TransientInput(name); - transientinput->AddTimeInput(new PentaInput(name,values,P1Enum),time); - transientinput->Configure(parameters); - } - this->inputs->AddInput(transientinput); - break; - - default: - _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet"); - } - -} -/*}}}*/ #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 20810) @@ -41,6 +41,6 @@ /*Object virtual functions definitions: {{{*/ Object *copy(); + void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction); int ObjectEnum(); - void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction); /*}}}*/ /*Penta routines:{{{*/ @@ -48,6 +48,6 @@ void AddInput(int input_enum, IssmDouble* values, int interpolation_enum); void AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part); + void CalvingRateLevermann(); IssmDouble CharacteristicLength(void){_error_("not implemented yet");}; - void CalvingRateLevermann(); void ComputeBasalStress(Vector* sigma_b); void ComputeDeviatoricStressTensor(); @@ -63,4 +63,5 @@ IssmDouble FloatingArea(void); void FSContactMigration(Vector* vertexgrounded,Vector* vertexfloating); + IssmDouble GetArea3D(void){_error_("not implemented yet!");}; void GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints); Element* GetBasalElement(void); @@ -172,9 +173,8 @@ void ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum); - IssmDouble GetArea3D(void){_error_("not implemented yet!");}; #ifdef _HAVE_DAKOTA_ + void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type); void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type); - void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type); #endif @@ -183,8 +183,8 @@ #endif #ifdef _HAVE_SEALEVELRISE_ + IssmDouble OceanArea(void){_error_("not implemented yet!");}; + IssmDouble OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");}; void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; - IssmDouble OceanArea(void){_error_("not implemented yet!");}; - IssmDouble OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");}; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 20810) @@ -190,4 +190,5 @@ IssmDouble y1,y2,y3,y4,y5,y6; IssmDouble z1,z2,z3,z4,z5,z6; + IssmDouble j_const_reciprocal; // SQRT3/12.0 /*Cast gauss to GaussPenta*/ @@ -224,15 +225,22 @@ z6=xyz_list[3*5+2]; + /* DEBUG + * REMOVE BEFORE SUBMITTING + * ---------- + * CHANGE SQRT3/12 TO MULTIPLY BY RECIPROCAL + */ + j_const_reciprocal=SQRT3/12; + J[NDOF3*0+0] = 0.25*(x1-x2-x4+x5)*zi+0.25*(-x1+x2-x4+x5); - J[NDOF3*1+0] = SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*zi+SQRT3/12.0*(-x1-x2+2*x3-x4-x5+2*x6); - J[NDOF3*2+0] = SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*eta+1/4*(x1-x2-x4+x5)*xi +0.25*(-x1+x5-x2+x4); + J[NDOF3*1+0] = j_const_reciprocal*(x1+x2-2*x3-x4-x5+2*x6)*zi+j_const_reciprocal*(-x1-x2+2*x3-x4-x5+2*x6); + J[NDOF3*2+0] = j_const_reciprocal*(x1+x2-2*x3-x4-x5+2*x6)*eta+0.25*(x1-x2-x4+x5)*xi +0.25*(-x1+x5-x2+x4); J[NDOF3*0+1] = 0.25*(y1-y2-y4+y5)*zi+0.25*(-y1+y2-y4+y5); - J[NDOF3*1+1] = SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*zi+SQRT3/12.0*(-y1-y2+2*y3-y4-y5+2*y6); - J[NDOF3*2+1] = SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*eta+0.25*(y1-y2-y4+y5)*xi+0.25*(y4-y1+y5-y2); + J[NDOF3*1+1] = j_const_reciprocal*(y1+y2-2*y3-y4-y5+2*y6)*zi+j_const_reciprocal*(-y1-y2+2*y3-y4-y5+2*y6); + J[NDOF3*2+1] = j_const_reciprocal*(y1+y2-2*y3-y4-y5+2*y6)*eta+0.25*(y1-y2-y4+y5)*xi+0.25*(y4-y1+y5-y2); J[NDOF3*0+2] = 0.25*(z1-z2-z4+z5)*zi+0.25*(-z1+z2-z4+z5); - J[NDOF3*1+2] = SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*zi+SQRT3/12.0*(-z1-z2+2*z3-z4-z5+2*z6); - J[NDOF3*2+2] = SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*eta+0.25*(z1-z2-z4+z5)*xi+0.25*(-z1+z5-z2+z4); + J[NDOF3*1+2] = j_const_reciprocal*(z1+z2-2*z3-z4-z5+2*z6)*zi+j_const_reciprocal*(-z1-z2+2*z3-z4-z5+2*z6); + J[NDOF3*2+2] = j_const_reciprocal*(z1+z2-2*z3-z4-z5+2*z6)*eta+0.25*(z1-z2-z4+z5)*xi+0.25*(-z1+z5-z2+z4); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h (revision 20810) @@ -28,4 +28,5 @@ void GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss); void GetTriaJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss); + void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*do nothing */}; int NumberofNodes(int finiteelement); int PressureInterpolation(int fe_stokes); @@ -33,5 +34,4 @@ int TensorInterpolation(int fe_stokes); int VelocityInterpolation(int fe_stokes); - void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*do nothing */}; }; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 20810) @@ -36,7 +36,7 @@ /*}}}*/ /*Object virtual functions definitions:{{{ */ - int ObjectEnum(); Object *copy(); void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction); + int ObjectEnum(); /*}}}*/ /*Element virtual functions definitions: {{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/SegRef.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/SegRef.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Elements/SegRef.h (revision 20810) @@ -24,6 +24,6 @@ void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussSeg* gauss,int finiteelement); void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussSeg* gauss,int finiteelement); + void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*do nothing */}; int NumberofNodes(int finiteelement); - void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*do nothing */}; }; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 20810) @@ -36,7 +36,7 @@ /*}}}*/ /*Object virtual functions definitions:{{{ */ - int ObjectEnum(); Object *copy(); void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction); + int ObjectEnum(); /*}}}*/ /*Element virtual functions definitions: {{{*/ @@ -62,4 +62,5 @@ IssmDouble FloatingArea(void){_error_("not implemented yet");}; void FSContactMigration(Vector* vertexgrounded,Vector* vertexfloating){_error_("not implemented yet");}; + IssmDouble GetArea3D(void){_error_("not implemented yet!");}; Element* GetBasalElement(void){_error_("not implemented yet");}; int GetElementType(void); @@ -169,5 +170,4 @@ void ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum); - IssmDouble GetArea3D(void){_error_("not implemented yet!");}; #ifdef _HAVE_GIA_ Index: /issm/trunk-jpl/src/c/classes/Elements/TetraRef.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/TetraRef.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Elements/TetraRef.h (revision 20810) @@ -25,9 +25,9 @@ void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTetra* gauss,int finiteelement); void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTetra* gauss,int finiteelement); + void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*do nothing */}; int NumberofNodes(int finiteelement); int PressureInterpolation(int fe_stokes); int TensorInterpolation(int fe_stokes); int VelocityInterpolation(int fe_stokes); - void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*do nothing */}; }; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 20810) @@ -208,59 +208,4 @@ } /*}}}*/ -void Tria::CalvingRateLevermann(){/*{{{*/ - - IssmDouble xyz_list[NUMVERTICES][3]; - GaussTria* gauss=NULL; - IssmDouble vx,vy,vel; - IssmDouble strainparallel; - IssmDouble propcoeff; - IssmDouble strainperpendicular; - IssmDouble calvingratex[NUMVERTICES]; - IssmDouble calvingratey[NUMVERTICES]; - IssmDouble calvingrate[NUMVERTICES]; - - - /* Get node coordinates and dof list: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - /*Retrieve all inputs and parameters we will need*/ - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* strainparallel_input=inputs->GetInput(StrainRateparallelEnum); _assert_(strainparallel_input); - Input* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum); _assert_(strainperpendicular_input); - Input* levermanncoeff_input=inputs->GetInput(CalvinglevermannCoeffEnum); _assert_(levermanncoeff_input); - - /* Start looping on the number of vertices: */ - gauss=new GaussTria(); - for (int iv=0;ivGaussVertex(iv); - - /* Get the value we need*/ - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - vel=vx*vx+vy*vy; - strainparallel_input->GetInputValue(&strainparallel,gauss); - strainperpendicular_input->GetInputValue(&strainperpendicular,gauss); - levermanncoeff_input->GetInputValue(&propcoeff,gauss); - - /*Calving rate proportionnal to the positive product of the strain rate along the ice flow direction and the strain rate perpendicular to the ice flow */ - calvingrate[iv]=propcoeff*strainparallel*strainperpendicular; - if(calvingrate[iv]<0){ - calvingrate[iv]=0; - } - calvingratex[iv]=calvingrate[iv]*vx/(sqrt(vel)+1.e-14); - calvingratey[iv]=calvingrate[iv]*vy/(sqrt(vel)+1.e-14); - } - - /*Add input*/ - this->inputs->AddInput(new TriaInput(CalvingratexEnum,&calvingratex[0],P1Enum)); - this->inputs->AddInput(new TriaInput(CalvingrateyEnum,&calvingratey[0],P1Enum)); - this->inputs->AddInput(new TriaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum)); - - /*Clean up and return*/ - delete gauss; - -} -/*}}}*/ void Tria::CalvingRateDev(){/*{{{*/ @@ -338,122 +283,57 @@ } /*}}}*/ -void Tria::WriteLevelsetSegment(DataSet* segments){/*{{{*/ - - if(!this->IsZeroLevelset(MaskIceLevelsetEnum)) return; - - IssmDouble* xyz_list_zero = NULL; +void Tria::CalvingRateLevermann(){/*{{{*/ + IssmDouble xyz_list[NUMVERTICES][3]; - ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES); - this->ZeroLevelsetCoordinates(&xyz_list_zero,&xyz_list[0][0], MaskIceLevelsetEnum); - if(xyz_list_zero){ - IssmDouble x[2]; - IssmDouble y[2]; - x[0] = xyz_list_zero[0*3 + 0]; x[1] = xyz_list_zero[1*3 + 0]; - y[0] = xyz_list_zero[0*3 + 1]; y[1] = xyz_list_zero[1*3 + 1]; - segments->AddObject(new Contour(segments->Size()+1,2,&x[0],&y[0],false)); - } - xDelete(xyz_list_zero); - -// IssmDouble ls[NUMVERTICES]; -// IssmDouble xyz_list[NUMVERTICES][3]; -// -// if(IsIceInElement()){ -// -// /*Retrieve all inputs and parameters*/ -// GetInputListOnVertices(&ls[0],levelset_enum); -// -// /*If the level set is awlays <0, there is no ice front here*/ -// bool iszerols= false; -// if(IsIceInElement()){ -// if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]*ls[2]==0. && ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]<=0.)){ -// iszerols = true; -// } -// } -// -// if(iszerols){ -// /*OK we have one segment!*/ -// ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); -// int count = 0; -// IssmDouble x[2]; -// IssmDouble y[2]; -// -// for(int i=0;i=0){ -// -// } -// -// } -// Contour* segment = new Contour(segment->Size()+1,2,x,y,false); -// } -// -// } -// -// _error_("STOP"); -} -/*}}}*/ -void Tria::ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments){/*{{{*/ - - /*Intermediaries*/ - IssmDouble d,xn,yn; - - /*Get current levelset and vertex coordinates*/ - IssmDouble ls[NUMVERTICES]; - IssmDouble xyz_list[NUMVERTICES][3]; + GaussTria* gauss=NULL; + IssmDouble vx,vy,vel; + IssmDouble strainparallel; + IssmDouble propcoeff; + IssmDouble strainperpendicular; + IssmDouble calvingratex[NUMVERTICES]; + IssmDouble calvingratey[NUMVERTICES]; + IssmDouble calvingrate[NUMVERTICES]; + + + /* Get node coordinates and dof list: */ ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum); - InputDuplicate(MaskIceLevelsetEnum,PressureEnum); - - /*Get distance from list of segments and reset ls*/ - for(int j=0;j 1.0){ - // Beyond the 'w' end of the segment - d = (x-segments[4*i+2])*(x-segments[4*i+2])+(y-segments[4*i+3])*(y-segments[4*i+3]); - } - else{ - // Projection falls on the segment - xn = segments[4*i+0] + t * (segments[4*i+2] - segments[4*i+0]); - yn = segments[4*i+1] + t * (segments[4*i+3] - segments[4*i+1]); - d = (x-xn)*(x-xn)+(y-yn)*(y-yn); - } - - if(d10000) dmin=10000; - if(ls[j]>0){ - ls[j] = dmin; - } - else{ - ls[j] = - dmin; - } - } - - /*Update Levelset*/ - this->inputs->AddInput(new TriaInput(MaskIceLevelsetEnum,&ls[0],P1Enum)); + + /*Retrieve all inputs and parameters we will need*/ + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* strainparallel_input=inputs->GetInput(StrainRateparallelEnum); _assert_(strainparallel_input); + Input* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum); _assert_(strainperpendicular_input); + Input* levermanncoeff_input=inputs->GetInput(CalvinglevermannCoeffEnum); _assert_(levermanncoeff_input); + + /* Start looping on the number of vertices: */ + gauss=new GaussTria(); + for (int iv=0;ivGaussVertex(iv); + + /* Get the value we need*/ + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + vel=vx*vx+vy*vy; + strainparallel_input->GetInputValue(&strainparallel,gauss); + strainperpendicular_input->GetInputValue(&strainperpendicular,gauss); + levermanncoeff_input->GetInputValue(&propcoeff,gauss); + + /*Calving rate proportionnal to the positive product of the strain rate along the ice flow direction and the strain rate perpendicular to the ice flow */ + calvingrate[iv]=propcoeff*strainparallel*strainperpendicular; + if(calvingrate[iv]<0){ + calvingrate[iv]=0; + } + calvingratex[iv]=calvingrate[iv]*vx/(sqrt(vel)+1.e-14); + calvingratey[iv]=calvingrate[iv]*vy/(sqrt(vel)+1.e-14); + } + + /*Add input*/ + this->inputs->AddInput(new TriaInput(CalvingratexEnum,&calvingratex[0],P1Enum)); + this->inputs->AddInput(new TriaInput(CalvingrateyEnum,&calvingratey[0],P1Enum)); + this->inputs->AddInput(new TriaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum)); + + /*Clean up and return*/ + delete gauss; + } /*}}}*/ @@ -2744,4 +2624,68 @@ } /*}}}*/ +void Tria::ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments){/*{{{*/ + + /*Intermediaries*/ + IssmDouble d,xn,yn; + + /*Get current levelset and vertex coordinates*/ + IssmDouble ls[NUMVERTICES]; + IssmDouble xyz_list[NUMVERTICES][3]; + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum); + InputDuplicate(MaskIceLevelsetEnum,PressureEnum); + + /*Get distance from list of segments and reset ls*/ + for(int j=0;j 1.0){ + // Beyond the 'w' end of the segment + d = (x-segments[4*i+2])*(x-segments[4*i+2])+(y-segments[4*i+3])*(y-segments[4*i+3]); + } + else{ + // Projection falls on the segment + xn = segments[4*i+0] + t * (segments[4*i+2] - segments[4*i+0]); + yn = segments[4*i+1] + t * (segments[4*i+3] - segments[4*i+1]); + d = (x-xn)*(x-xn)+(y-yn)*(y-yn); + } + + if(d10000) dmin=10000; + if(ls[j]>0){ + ls[j] = dmin; + } + else{ + ls[j] = - dmin; + } + } + + /*Update Levelset*/ + this->inputs->AddInput(new TriaInput(MaskIceLevelsetEnum,&ls[0],P1Enum)); +} +/*}}}*/ void Tria::SetClone(int* minranks){/*{{{*/ @@ -3361,4 +3305,60 @@ } /*}}}*/ +void Tria::WriteLevelsetSegment(DataSet* segments){/*{{{*/ + + if(!this->IsZeroLevelset(MaskIceLevelsetEnum)) return; + + IssmDouble* xyz_list_zero = NULL; + IssmDouble xyz_list[NUMVERTICES][3]; + ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES); + this->ZeroLevelsetCoordinates(&xyz_list_zero,&xyz_list[0][0], MaskIceLevelsetEnum); + if(xyz_list_zero){ + IssmDouble x[2]; + IssmDouble y[2]; + x[0] = xyz_list_zero[0*3 + 0]; x[1] = xyz_list_zero[1*3 + 0]; + y[0] = xyz_list_zero[0*3 + 1]; y[1] = xyz_list_zero[1*3 + 1]; + segments->AddObject(new Contour(segments->Size()+1,2,&x[0],&y[0],false)); + } + xDelete(xyz_list_zero); + +// IssmDouble ls[NUMVERTICES]; +// IssmDouble xyz_list[NUMVERTICES][3]; +// +// if(IsIceInElement()){ +// +// /*Retrieve all inputs and parameters*/ +// GetInputListOnVertices(&ls[0],levelset_enum); +// +// /*If the level set is awlays <0, there is no ice front here*/ +// bool iszerols= false; +// if(IsIceInElement()){ +// if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]*ls[2]==0. && ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]<=0.)){ +// iszerols = true; +// } +// } +// +// if(iszerols){ +// /*OK we have one segment!*/ +// ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); +// int count = 0; +// IssmDouble x[2]; +// IssmDouble y[2]; +// +// for(int i=0;i=0){ +// +// } +// +// } +// Contour* segment = new Contour(segment->Size()+1,2,x,y,false); +// } +// +// } +// +// _error_("STOP"); +} +/*}}}*/ void Tria::ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){/*{{{*/ /* Return coordinates where levelset intersects element edges. @@ -3541,4 +3541,30 @@ #ifdef _HAVE_SEALEVELRISE_ +IssmDouble Tria::OceanArea(void){ /*{{{*/ + + if(IsWaterInElement()) return GetArea3D(); + else return 0; + +} +/*}}}*/ +IssmDouble Tria::OceanAverage(IssmDouble* Sg){ /*{{{*/ + + if(IsWaterInElement()){ + + IssmDouble area; + + /*Compute area of element:*/ + area=GetArea3D(); + + /*Average Sg over vertices:*/ + IssmDouble Sg_avg=0; for(int i=0;ivertices[i]->Sid()]/NUMVERTICES; + + /*return: */ + return area*Sg_avg; + } + else return 0; + +} +/*}}}*/ void Tria::SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/ @@ -3827,30 +3853,4 @@ return; -} -/*}}}*/ -IssmDouble Tria::OceanAverage(IssmDouble* Sg){ /*{{{*/ - - if(IsWaterInElement()){ - - IssmDouble area; - - /*Compute area of element:*/ - area=GetArea3D(); - - /*Average Sg over vertices:*/ - IssmDouble Sg_avg=0; for(int i=0;ivertices[i]->Sid()]/NUMVERTICES; - - /*return: */ - return area*Sg_avg; - } - else return 0; - -} -/*}}}*/ -IssmDouble Tria::OceanArea(void){ /*{{{*/ - - if(IsWaterInElement()) return GetArea3D(); - else return 0; - } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 20810) @@ -38,7 +38,7 @@ /*}}}*/ /*Object virtual functions definitions:{{{ */ - int ObjectEnum(); Object *copy(); void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction); + int ObjectEnum(); /*}}}*/ /*Update virtual functions resolution: {{{*/ @@ -52,8 +52,6 @@ /*Element virtual functions definitions: {{{*/ void AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part); + void CalvingRateDev(); void CalvingRateLevermann(); - void CalvingRateDev(); - void WriteLevelsetSegment(DataSet* segments); - void ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments); IssmDouble CharacteristicLength(void); void ComputeBasalStress(Vector* sigma_b); @@ -119,4 +117,5 @@ void ResetFSBasalBoundaryCondition(void); void ResetHooks(); + void ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments); void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index); void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters); @@ -139,4 +138,5 @@ int VertexConnectivity(int vertexindex); void VerticalSegmentIndices(int** pindices,int* pnumseg){_error_("not implemented yet");}; + void WriteLevelsetSegment(DataSet* segments); void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum); @@ -145,8 +145,8 @@ #endif #ifdef _HAVE_SEALEVELRISE_ + IssmDouble OceanArea(void); + IssmDouble OceanAverage(IssmDouble* Sg); void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea); void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea); - IssmDouble OceanAverage(IssmDouble* Sg); - IssmDouble OceanArea(void); #endif /*}}}*/ @@ -158,7 +158,7 @@ IssmDouble GetAreaIce(void); void GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints); - void GetLevelsetIntersection(int** pindices, int* pnumiceverts, IssmDouble* fraction, int levelset_enum, IssmDouble level); int GetElementType(void); void GetInputValue(IssmDouble* pvalue,Node* node,int enumtype); + void GetLevelsetIntersection(int** pindices, int* pnumiceverts, IssmDouble* fraction, int levelset_enum, IssmDouble level); void GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype); Node* GetNode(int node_number); Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 20810) @@ -212,4 +212,107 @@ } +} +/*}}}*/ +void TriaRef::GetSegmentBFlux(IssmDouble* B,Gauss* gauss, int index1,int index2,int finiteelement){/*{{{*/ + /*Compute B matrix. B=[phi1 phi2 -phi3 -phi4] + * + * and phi1=phi3 phi2=phi4 + * + * We assume B has been allocated already, of size: 1x4 + */ + + /*Fetch number of nodes for this finite element*/ + int numnodes = this->NumberofNodes(finiteelement); + + /*Get nodal functions*/ + IssmDouble* basis=xNew(numnodes); + GetNodalFunctions(basis,gauss,finiteelement); + + /*Build B for this segment*/ + B[0] = +basis[index1]; + B[1] = +basis[index2]; + B[2] = -basis[index1]; + B[3] = -basis[index2]; + + /*Clean-up*/ + xDelete(basis); +} +/*}}}*/ +void TriaRef::GetSegmentBprimeFlux(IssmDouble* Bprime,Gauss* gauss, int index1,int index2,int finiteelement){/*{{{*/ + /*Compute Bprime matrix. Bprime=[phi1 phi2 phi3 phi4] + * + * and phi1=phi3 phi2=phi4 + * + * We assume Bprime has been allocated already, of size: 1x4 + */ + + /*Fetch number of nodes for this finite element*/ + int numnodes = this->NumberofNodes(finiteelement); + + /*Get nodal functions*/ + IssmDouble* basis=xNew(numnodes); + GetNodalFunctions(basis,gauss,finiteelement); + + /*Build B'*/ + Bprime[0] = basis[index1]; + Bprime[1] = basis[index2]; + Bprime[2] = basis[index1]; + Bprime[3] = basis[index2]; + + /*Clean-up*/ + xDelete(basis); +} +/*}}}*/ +void TriaRef::GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + /*The Jacobian determinant is constant over the element, discard the gaussian points. + * J is assumed to have been allocated*/ + + IssmDouble x1 = xyz_list[3*0+0]; + IssmDouble y1 = xyz_list[3*0+1]; + IssmDouble x2 = xyz_list[3*1+0]; + IssmDouble y2 = xyz_list[3*1+1]; + + *Jdet = .5*sqrt(pow(x2-x1,2) + pow(y2-y1,2)); + if(*Jdet<0) _error_("negative jacobian determinant!"); + +} +/*}}}*/ +void TriaRef::GetSegmentNodalFunctions(IssmDouble* basis,Gauss* gauss,int index1,int index2,int finiteelement){/*{{{*/ + /*This routine returns the values of the nodal functions at the gaussian point.*/ + + _assert_(index1>=0 && index1<3); + _assert_(index2>=0 && index2<3); + + /*Fetch number of nodes for this finite element*/ + int numnodes = this->NumberofNodes(finiteelement); + + /*Get nodal functions*/ + IssmDouble* triabasis=xNew(numnodes); + GetNodalFunctions(triabasis,gauss,finiteelement); + + switch(finiteelement){ + case P1Enum: case P1DGEnum: + basis[0]=triabasis[index1]; + basis[1]=triabasis[index2]; + xDelete(triabasis); + return; + case P1bubbleEnum: case P1bubblecondensedEnum: + basis[0]=triabasis[index1]; + basis[1]=triabasis[index2]; + xDelete(triabasis); + return; + case P2Enum: + _assert_(index2(triabasis); + return; + default: + _error_("Element type "<(triabasis); } /*}}}*/ @@ -302,107 +405,4 @@ } -} -/*}}}*/ -void TriaRef::GetSegmentBFlux(IssmDouble* B,Gauss* gauss, int index1,int index2,int finiteelement){/*{{{*/ - /*Compute B matrix. B=[phi1 phi2 -phi3 -phi4] - * - * and phi1=phi3 phi2=phi4 - * - * We assume B has been allocated already, of size: 1x4 - */ - - /*Fetch number of nodes for this finite element*/ - int numnodes = this->NumberofNodes(finiteelement); - - /*Get nodal functions*/ - IssmDouble* basis=xNew(numnodes); - GetNodalFunctions(basis,gauss,finiteelement); - - /*Build B for this segment*/ - B[0] = +basis[index1]; - B[1] = +basis[index2]; - B[2] = -basis[index1]; - B[3] = -basis[index2]; - - /*Clean-up*/ - xDelete(basis); -} -/*}}}*/ -void TriaRef::GetSegmentBprimeFlux(IssmDouble* Bprime,Gauss* gauss, int index1,int index2,int finiteelement){/*{{{*/ - /*Compute Bprime matrix. Bprime=[phi1 phi2 phi3 phi4] - * - * and phi1=phi3 phi2=phi4 - * - * We assume Bprime has been allocated already, of size: 1x4 - */ - - /*Fetch number of nodes for this finite element*/ - int numnodes = this->NumberofNodes(finiteelement); - - /*Get nodal functions*/ - IssmDouble* basis=xNew(numnodes); - GetNodalFunctions(basis,gauss,finiteelement); - - /*Build B'*/ - Bprime[0] = basis[index1]; - Bprime[1] = basis[index2]; - Bprime[2] = basis[index1]; - Bprime[3] = basis[index2]; - - /*Clean-up*/ - xDelete(basis); -} -/*}}}*/ -void TriaRef::GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*The Jacobian determinant is constant over the element, discard the gaussian points. - * J is assumed to have been allocated*/ - - IssmDouble x1 = xyz_list[3*0+0]; - IssmDouble y1 = xyz_list[3*0+1]; - IssmDouble x2 = xyz_list[3*1+0]; - IssmDouble y2 = xyz_list[3*1+1]; - - *Jdet = .5*sqrt(pow(x2-x1,2) + pow(y2-y1,2)); - if(*Jdet<0) _error_("negative jacobian determinant!"); - -} -/*}}}*/ -void TriaRef::GetSegmentNodalFunctions(IssmDouble* basis,Gauss* gauss,int index1,int index2,int finiteelement){/*{{{*/ - /*This routine returns the values of the nodal functions at the gaussian point.*/ - - _assert_(index1>=0 && index1<3); - _assert_(index2>=0 && index2<3); - - /*Fetch number of nodes for this finite element*/ - int numnodes = this->NumberofNodes(finiteelement); - - /*Get nodal functions*/ - IssmDouble* triabasis=xNew(numnodes); - GetNodalFunctions(triabasis,gauss,finiteelement); - - switch(finiteelement){ - case P1Enum: case P1DGEnum: - basis[0]=triabasis[index1]; - basis[1]=triabasis[index2]; - xDelete(triabasis); - return; - case P1bubbleEnum: case P1bubblecondensedEnum: - basis[0]=triabasis[index1]; - basis[1]=triabasis[index2]; - xDelete(triabasis); - return; - case P2Enum: - _assert_(index2(triabasis); - return; - default: - _error_("Element type "<(triabasis); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h (revision 20810) @@ -28,4 +28,5 @@ void GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss); void GetSegmentNodalFunctions(IssmDouble* basis,Gauss* gauss, int index1,int index2,int finiteelement); + void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*do nothing */}; void NodeOnEdgeIndices(int* pnumindices,int** pindices,int index,int finiteelement); int NumberofNodes(int finiteelement); @@ -33,5 +34,4 @@ int TensorInterpolation(int fe_stokes); int VelocityInterpolation(int fe_stokes); - void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*do nothing */}; }; Index: /issm/trunk-jpl/src/c/classes/ExternalResults/ExternalResult.h =================================================================== --- /issm/trunk-jpl/src/c/classes/ExternalResults/ExternalResult.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/ExternalResults/ExternalResult.h (revision 20810) @@ -24,9 +24,9 @@ virtual ~ExternalResult(){}; - virtual void WriteData(FILE* fid,bool io_gather)=0; + virtual int GetResultEnum(void)=0; virtual char* GetResultName(void)=0; virtual int GetStep(void)=0; virtual double GetValue(void)=0; - virtual int GetResultEnum(void)=0; + virtual void WriteData(FILE* fid,bool io_gather)=0; }; #endif Index: /issm/trunk-jpl/src/c/classes/ExternalResults/GenericExternalResult.h =================================================================== --- /issm/trunk-jpl/src/c/classes/ExternalResults/GenericExternalResult.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/ExternalResults/GenericExternalResult.h (revision 20810) @@ -132,4 +132,7 @@ /*Object virtual functions definitions:*/ + Object* copy(void) { /*{{{*/ + return new GenericExternalResult(this->id,this->result_name,this->value,this->step,this->time); + } /*}}}*/ void Echo(void){ /*{{{*/ this->DeepEcho(); @@ -146,7 +149,4 @@ _error_("template ObjectEnum not implemented for this ResultType\n"); } /*}}}*/ - Object* copy(void) { /*{{{*/ - return new GenericExternalResult(this->id,this->result_name,this->value,this->step,this->time); - } /*}}}*/ void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){/*{{{*/ _error_("not implemented yet!"); Index: /issm/trunk-jpl/src/c/classes/ExternalResults/Results.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/ExternalResults/Results.cpp (revision 20809) +++ /issm/trunk-jpl/src/c/classes/ExternalResults/Results.cpp (revision 20810) @@ -31,20 +31,4 @@ /*Object management*/ -void Results::Write(Parameters* parameters){/*{{{*/ - - FILE *fid = NULL; - bool io_gather; - - /*Recover file descriptor: */ - parameters->FindParam(&fid,OutputFilePointerEnum); - parameters->FindParam(&io_gather,SettingsIoGatherEnum); - - for(int i=0;iSize();i++){ - ExternalResult* result=xDynamicCast(this->GetObjectByOffset(i)); - result->WriteData(fid,io_gather); - } - -} -/*}}}*/ int Results::AddResult(ExternalResult* in_result){/*{{{*/ @@ -106,2 +90,18 @@ } /*}}}*/ +void Results::Write(Parameters* parameters){/*{{{*/ + + FILE *fid = NULL; + bool io_gather; + + /*Recover file descriptor: */ + parameters->FindParam(&fid,OutputFilePointerEnum); + parameters->FindParam(&io_gather,SettingsIoGatherEnum); + + for(int i=0;iSize();i++){ + ExternalResult* result=xDynamicCast(this->GetObjectByOffset(i)); + result->WriteData(fid,io_gather); + } + +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/ExternalResults/Results.h =================================================================== --- /issm/trunk-jpl/src/c/classes/ExternalResults/Results.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/ExternalResults/Results.h (revision 20810) @@ -22,6 +22,6 @@ /*Mehthos*/ int AddResult(ExternalResult* result); + int DeleteResult(int result_enum,int result_step); ExternalResult* FindResult(int result_enum); - int DeleteResult(int result_enum,int result_step); void Write(Parameters* parameters); }; Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 20809) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 20810) @@ -143,90 +143,45 @@ /*Object management*/ -void FemModel::Echo(void){/*{{{*/ - - _printf_("FemModel echo: \n"); - _printf_(" number of fem models: " << nummodels << "\n"); - _printf_(" analysis_type_list: \n"); - for(int i=0;isolution_type = in_solution_type; - this->analysis_counter = nummodels-1; //point to last analysis_type carried out. - this->results = new Results(); //not initialized by CreateDataSets - /*Open input file on cpu 0 and create IoModel */ - if(my_rank==0) IOMODEL = pfopen0(inputfilename ,"rb"); - IoModel* iomodel = new IoModel(IOMODEL,in_solution_type,trace,X); - - /*Figure out what analyses are activated for this solution*/ - SolutionAnalysesList(&this->analysis_type_list,&this->nummodels,iomodel,this->solution_type); - - /*Open toolkits file: */ - toolkitsoptionsfid=pfopen(toolkitsfilename,"r"); - - /*create datasets for all analyses*/ - ModelProcessorx(&this->elements,&this->nodes,&this->vertices,&this->materials,&this->constraints,&this->loads,&this->parameters,iomodel,toolkitsoptionsfid,rootpath,this->solution_type,this->nummodels,this->analysis_type_list); - - - /*do the post-processing of the datasets to get an FemModel that can actually run analyses: */ - for(i=0;iSetCurrentConfiguration(analysis_type); - - if(VerboseMProcessor()) _printf0_(" configuring element and loads\n"); - ConfigureObjectsx(elements, loads, nodes, vertices, materials,parameters); - - if(i==0){ - if(VerboseMProcessor()) _printf0_(" creating vertex PIDs\n"); - VerticesDofx(vertices,parameters); - - if(VerboseMProcessor()) _printf0_(" detecting active vertices\n"); - GetMaskOfIceVerticesLSMx(this); - } - - if(VerboseMProcessor()) _printf0_(" resolving node constraints\n"); - SpcNodesx(nodes,constraints,parameters,analysis_type); - - if(VerboseMProcessor()) _printf0_(" creating nodal degrees of freedom\n"); - NodesDofx(nodes,parameters,analysis_type); - } - - /*Close input file and toolkits file descriptors: */ - if(my_rank==0) pfclose(IOMODEL,inputfilename); - pfclose(toolkitsoptionsfid,toolkitsfilename); - - /*Open output file once for all and add output file descriptor to parameters*/ - output_fid=pfopen(outputfilename,"wb"); - this->parameters->SetParam(output_fid,OutputFilePointerEnum); +void FemModel::CheckPoint(void){/*{{{*/ + + FILE* restartfid=NULL; + char* restartfilename = NULL; + int femmodel_size; + char* femmodel_buffer=NULL; + char* femmodel_buffer_ini=NULL; + + /*First, recover the name of the restart file: */ + parameters->FindParam(&restartfilename,RestartFileNameEnum); - /*Now save all of these file names into parameters, you never know when you might need them: */ - this->parameters->AddObject(new StringParam(ToolkitsFileNameEnum,toolkitsfilename)); - this->parameters->AddObject(new StringParam(RootPathEnum,rootpath)); - this->parameters->AddObject(new StringParam(InputFileNameEnum,inputfilename)); - this->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename)); - this->parameters->AddObject(new StringParam(LockFileNameEnum,lockfilename)); - this->parameters->AddObject(new StringParam(RestartFileNameEnum,restartfilename)); - - /*Clean up*/ - delete iomodel; -}/*}}}*/ + /*Open file for writing: */ + restartfid=pfopen(restartfilename,"wb"); + + /*Initialize: */ + femmodel_size=0; + + /*Create buffer to hold marshalled femmodel: */ + this->Marshall(NULL,&femmodel_size,MARSHALLING_SIZE); + femmodel_buffer=xNew(femmodel_size); + /*Keep track of initial position of femmodel_buffer: */ + femmodel_buffer_ini=femmodel_buffer; + + /*Marshall:*/ + this->Marshall(&femmodel_buffer,NULL,MARSHALLING_FORWARD); + + /*Reset position of buffer: */ + femmodel_buffer=femmodel_buffer_ini; + + /*write buffer: */ + fwrite(femmodel_buffer,femmodel_size,sizeof(char),restartfid); + + /*Done, close file :*/ + pfclose(restartfid,restartfilename); + + /*Free ressources: */ + xDelete(femmodel_buffer); + xDelete(restartfilename); + +} +/*}}}*/ void FemModel::CleanUp(void){/*{{{*/ @@ -275,4 +230,260 @@ } /*}}}*/ +FemModel* FemModel::copy(void){/*{{{*/ + + FemModel* output=NULL; + int i; + int analysis_type; + + output=new FemModel(*this); //Use default copy constructor. + + output->nummodels = this->nummodels; + output->solution_type = this->solution_type; + output->analysis_counter = this->analysis_counter; + + /*Now, deep copy arrays: */ + output->analysis_type_list=xNew(nummodels); + xMemCpy(output->analysis_type_list,this->analysis_type_list,this->nummodels); + + output->profiler=static_cast(this->profiler->copy()); + + output->loads=static_cast(this->loads->Copy()); + output->materials=static_cast(this->materials->Copy()); + output->parameters=static_cast(this->parameters->Copy()); + output->constraints=static_cast(this->constraints->Copy()); + output->results=static_cast(this->results->Copy()); + + output->nodes=static_cast(this->nodes->Copy()); + output->vertices=static_cast(this->vertices->Copy()); + output->elements=static_cast(this->elements->Copy()); + + /*reset hooks for elements, loads and nodes: */ + output->elements->ResetHooks(); + output->loads->ResetHooks(); + output->materials->ResetHooks(); + + /*do the post-processing of the datasets to get an FemModel that can actually run analyses: */ + for(i=0;ianalysis_type_list[i]; + output->SetCurrentConfiguration(analysis_type); + if(i==0) VerticesDofx(output->vertices,output->parameters); //only call once, we only have one set of vertices + SpcNodesx(output->nodes,output->constraints,output->parameters,analysis_type); + NodesDofx(output->nodes,output->parameters,analysis_type); + ConfigureObjectsx(output->elements,output->loads,output->nodes,output->vertices,output->materials,output->parameters); + } + + /*Reset current configuration: */ + analysis_type=output->analysis_type_list[analysis_counter]; + output->SetCurrentConfiguration(analysis_type); + + return output; +} +/*}}}*/ +void FemModel::Echo(void){/*{{{*/ + + _printf_("FemModel echo: \n"); + _printf_(" number of fem models: " << nummodels << "\n"); + _printf_(" analysis_type_list: \n"); + for(int i=0;isolution_type = in_solution_type; + this->analysis_counter = nummodels-1; //point to last analysis_type carried out. + this->results = new Results(); //not initialized by CreateDataSets + /*Open input file on cpu 0 and create IoModel */ + if(my_rank==0) IOMODEL = pfopen0(inputfilename ,"rb"); + IoModel* iomodel = new IoModel(IOMODEL,in_solution_type,trace,X); + + /*Figure out what analyses are activated for this solution*/ + SolutionAnalysesList(&this->analysis_type_list,&this->nummodels,iomodel,this->solution_type); + + /*Open toolkits file: */ + toolkitsoptionsfid=pfopen(toolkitsfilename,"r"); + + /*create datasets for all analyses*/ + ModelProcessorx(&this->elements,&this->nodes,&this->vertices,&this->materials,&this->constraints,&this->loads,&this->parameters,iomodel,toolkitsoptionsfid,rootpath,this->solution_type,this->nummodels,this->analysis_type_list); + + + /*do the post-processing of the datasets to get an FemModel that can actually run analyses: */ + for(i=0;iSetCurrentConfiguration(analysis_type); + + if(VerboseMProcessor()) _printf0_(" configuring element and loads\n"); + ConfigureObjectsx(elements, loads, nodes, vertices, materials,parameters); + + if(i==0){ + if(VerboseMProcessor()) _printf0_(" creating vertex PIDs\n"); + VerticesDofx(vertices,parameters); + + if(VerboseMProcessor()) _printf0_(" detecting active vertices\n"); + GetMaskOfIceVerticesLSMx(this); + } + + if(VerboseMProcessor()) _printf0_(" resolving node constraints\n"); + SpcNodesx(nodes,constraints,parameters,analysis_type); + + if(VerboseMProcessor()) _printf0_(" creating nodal degrees of freedom\n"); + NodesDofx(nodes,parameters,analysis_type); + } + + /*Close input file and toolkits file descriptors: */ + if(my_rank==0) pfclose(IOMODEL,inputfilename); + pfclose(toolkitsoptionsfid,toolkitsfilename); + + /*Open output file once for all and add output file descriptor to parameters*/ + output_fid=pfopen(outputfilename,"wb"); + this->parameters->SetParam(output_fid,OutputFilePointerEnum); + + /*Now save all of these file names into parameters, you never know when you might need them: */ + this->parameters->AddObject(new StringParam(ToolkitsFileNameEnum,toolkitsfilename)); + this->parameters->AddObject(new StringParam(RootPathEnum,rootpath)); + this->parameters->AddObject(new StringParam(InputFileNameEnum,inputfilename)); + this->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename)); + this->parameters->AddObject(new StringParam(LockFileNameEnum,lockfilename)); + this->parameters->AddObject(new StringParam(RestartFileNameEnum,restartfilename)); + + /*Clean up*/ + delete iomodel; +}/*}}}*/ +void FemModel::Marshall(char** pmarshalled_data, int* pmarshalled_data_size, int marshall_direction){ /*{{{*/ + + int i; + int analysis_type; + + if(marshall_direction==MARSHALLING_BACKWARD){ + delete this->profiler; + delete this->loads; + delete this->materials; + delete this->parameters; + delete this->constraints; + delete this->results; + delete this->nodes; + delete this->vertices; + delete this->elements; + xDelete(this->analysis_type_list); + + this->profiler = new Profiler(); + this->loads = new Loads(); + this->materials = new Materials(); + this->parameters = new Parameters(); + this->constraints = new Constraints(); + this->results = new Results(); + this->nodes = new Nodes(); + this->vertices = new Vertices(); + this->elements = new Elements(); + } + + MARSHALLING_ENUM(FemModelEnum); + + MARSHALLING(solution_type); + MARSHALLING(analysis_counter); + MARSHALLING(nummodels); + MARSHALLING_DYNAMIC(analysis_type_list,int,nummodels); + + this->profiler->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); + this->loads->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); + this->materials->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); + this->parameters->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); + this->constraints->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); + this->results->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); + this->nodes->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); + this->vertices->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); + this->elements->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); + + if(marshall_direction==MARSHALLING_BACKWARD){ + //reset hooks for elements, loads and nodes: + this->elements->ResetHooks(); + this->loads->ResetHooks(); + this->materials->ResetHooks(); + + //do the post-processing of the datasets to get an FemModel that can actually run analyses: + for(i=0;ianalysis_type_list[i]; + SetCurrentConfiguration(analysis_type); + if(i==0) VerticesDofx(this->vertices,this->parameters); //only call once, we only have one set of vertices + SpcNodesx(this->nodes,this->constraints,this->parameters,analysis_type); + NodesDofx(this->nodes,this->parameters,analysis_type); + ConfigureObjectsx(this->elements,this->loads,this->nodes,this->vertices,this->materials,this->parameters); + } + + //Reset current configuration: + analysis_type=this->analysis_type_list[analysis_counter]; + SetCurrentConfiguration(analysis_type); + } + +} +/*}}}*/ +void FemModel::Restart(){ /*{{{*/ + + FILE* restartfid=NULL; + FILE* output_fid=NULL; + char* restartfilename = NULL; + int femmodel_size=0; + int fread_return=0; + char* femmodel_buffer=NULL; + char* femmodel_buffer_ini=NULL; + + /*First, recover the name of the restart file: */ + parameters->FindParam(&restartfilename,RestartFileNameEnum); + + /*Now, figure out whether this file actually exists!: */ + restartfid=pfopen(restartfilename,"r",false); + + if(restartfid==NULL){ + xDelete(restartfilename); + return; //could not find the file, so no restart possible. + } + + /*Figure out size of buffer to be read: */ + fseek(restartfid, 0L, SEEK_END); + femmodel_size = ftell(restartfid); + fseek(restartfid, 0L, SEEK_SET); + + /*Allocate buffer: */ + femmodel_buffer=xNew(femmodel_size); + + /*Read buffer from file: */ + fread_return=fread(femmodel_buffer,femmodel_size,sizeof(char),restartfid); if(fread_return!=1)_error_("error reading the buffer from marshalled file!"); + femmodel_buffer_ini=femmodel_buffer; //keep track of the initial position, so as to free later. + + this->parameters->FindParam(&output_fid,OutputFilePointerEnum); + + /*Create new FemModel by demarshalling the buffer: */ + this->Marshall(&femmodel_buffer,NULL,MARSHALLING_BACKWARD); + + this->parameters->SetParam(output_fid,OutputFilePointerEnum); + + /*Reset position of buffer: */ + femmodel_buffer=femmodel_buffer_ini; + + /*Done, close file :*/ + pfclose(restartfid,restartfilename); + + /*Free ressources: */ + xDelete(restartfilename); + xDelete(femmodel_buffer); + +} +/*}}}*/ void FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){/*{{{*/ @@ -311,115 +522,4 @@ void FemModel::SetCurrentConfiguration(int configuration_type){/*{{{*/ this->SetCurrentConfiguration(configuration_type,configuration_type); -} -/*}}}*/ -void FemModel::Solve(void){/*{{{*/ - - /*profiling: */ - bool profiling = false; - IssmDouble solution_time; - IssmDouble solution_flops; - IssmDouble solution_memory; - - /*solution: */ - int solution_type; - void (*solutioncore)(FemModel*)=NULL; //core solution function pointer - - _printf0_("call computational core:\n"); - - /*Retrieve solution_type from parameters: */ - parameters->FindParam(&solution_type,SolutionTypeEnum); - - /*Figure out which solution core we are going to run with the current solution type: */ - WrapperCorePointerFromSolutionEnum(&solutioncore,this->parameters,solution_type); - - /*run solution core: */ - profiler->Tag(STARTCORE); - solutioncore(this); - profiler->Tag(FINISHCORE); - - /*run AD core if needed: */ - profiler->Tag(STARTADCORE); - ad_core(this); - profiler->Tag(FINISHADCORE); - - /*some profiling results for the core: */ - parameters->FindParam(&profiling,DebugProfilingEnum); - if(profiling){ - - solution_time=profiler->DeltaTime(STARTCORE,FINISHCORE); - solution_flops=profiler->DeltaFlops(STARTCORE,FINISHCORE); - solution_memory=profiler->Memory(FINISHCORE); - - _printf0_("Core solution elapsed time : " << solution_time << " Seconds\n"); - _printf0_("Core solution number of flops : " << solution_flops << " Flops\n"); - _printf0_("Core solution memory used : " << solution_memory << " Bytes\n"); - - /*Add to results: */ - results->AddObject(new GenericExternalResult(results->Size()+1, ProfilingSolutionTimeEnum, solution_time, 0, 0)); - results->AddObject(new GenericExternalResult(results->Size()+1, ProfilingCurrentMemEnum, solution_memory, 0, 0)); - results->AddObject(new GenericExternalResult(results->Size()+1, ProfilingCurrentFlopsEnum, solution_flops, 0, 0)); - - #ifdef _HAVE_ADOLC_ - solution_time=profiler->DeltaTime(STARTADCORE,FINISHADCORE); - solution_flops=profiler->DeltaFlops(STARTADCORE,FINISHADCORE); - solution_memory=profiler->Memory(FINISHADCORE); - - _printf0_("AD Solution elapsed time : " << solution_time << " Seconds\n"); - _printf0_("AD Solution number of flops : " << solution_flops << " Flops\n"); - _printf0_("AD Solution memory used : " << solution_memory << " Bytes\n"); - - #endif - - } -} -/*}}}*/ -FemModel* FemModel::copy(void){/*{{{*/ - - FemModel* output=NULL; - int i; - int analysis_type; - - output=new FemModel(*this); //Use default copy constructor. - - output->nummodels = this->nummodels; - output->solution_type = this->solution_type; - output->analysis_counter = this->analysis_counter; - - /*Now, deep copy arrays: */ - output->analysis_type_list=xNew(nummodels); - xMemCpy(output->analysis_type_list,this->analysis_type_list,this->nummodels); - - output->profiler=static_cast(this->profiler->copy()); - - output->loads=static_cast(this->loads->Copy()); - output->materials=static_cast(this->materials->Copy()); - output->parameters=static_cast(this->parameters->Copy()); - output->constraints=static_cast(this->constraints->Copy()); - output->results=static_cast(this->results->Copy()); - - output->nodes=static_cast(this->nodes->Copy()); - output->vertices=static_cast(this->vertices->Copy()); - output->elements=static_cast(this->elements->Copy()); - - /*reset hooks for elements, loads and nodes: */ - output->elements->ResetHooks(); - output->loads->ResetHooks(); - output->materials->ResetHooks(); - - /*do the post-processing of the datasets to get an FemModel that can actually run analyses: */ - for(i=0;ianalysis_type_list[i]; - output->SetCurrentConfiguration(analysis_type); - if(i==0) VerticesDofx(output->vertices,output->parameters); //only call once, we only have one set of vertices - SpcNodesx(output->nodes,output->constraints,output->parameters,analysis_type); - NodesDofx(output->nodes,output->parameters,analysis_type); - ConfigureObjectsx(output->elements,output->loads,output->nodes,output->vertices,output->materials,output->parameters); - } - - /*Reset current configuration: */ - analysis_type=output->analysis_type_list[analysis_counter]; - output->SetCurrentConfiguration(analysis_type); - - return output; } /*}}}*/ @@ -634,232 +734,208 @@ else xDelete(analyses); }/*}}}*/ -void FemModel::CheckPoint(void){/*{{{*/ - - FILE* restartfid=NULL; - char* restartfilename = NULL; - int femmodel_size; - char* femmodel_buffer=NULL; - char* femmodel_buffer_ini=NULL; - - /*First, recover the name of the restart file: */ - parameters->FindParam(&restartfilename,RestartFileNameEnum); - - /*Open file for writing: */ - restartfid=pfopen(restartfilename,"wb"); - - /*Initialize: */ - femmodel_size=0; - - /*Create buffer to hold marshalled femmodel: */ - this->Marshall(NULL,&femmodel_size,MARSHALLING_SIZE); - femmodel_buffer=xNew(femmodel_size); - /*Keep track of initial position of femmodel_buffer: */ - femmodel_buffer_ini=femmodel_buffer; - - /*Marshall:*/ - this->Marshall(&femmodel_buffer,NULL,MARSHALLING_FORWARD); - - /*Reset position of buffer: */ - femmodel_buffer=femmodel_buffer_ini; - - /*write buffer: */ - fwrite(femmodel_buffer,femmodel_size,sizeof(char),restartfid); - - /*Done, close file :*/ - pfclose(restartfid,restartfilename); - - /*Free ressources: */ - xDelete(femmodel_buffer); - xDelete(restartfilename); - -} -/*}}}*/ -void FemModel::Restart(){ /*{{{*/ - - FILE* restartfid=NULL; - FILE* output_fid=NULL; - char* restartfilename = NULL; - int femmodel_size=0; - int fread_return=0; - char* femmodel_buffer=NULL; - char* femmodel_buffer_ini=NULL; - - /*First, recover the name of the restart file: */ - parameters->FindParam(&restartfilename,RestartFileNameEnum); - - /*Now, figure out whether this file actually exists!: */ - restartfid=pfopen(restartfilename,"r",false); - - if(restartfid==NULL){ - xDelete(restartfilename); - return; //could not find the file, so no restart possible. - } - - /*Figure out size of buffer to be read: */ - fseek(restartfid, 0L, SEEK_END); - femmodel_size = ftell(restartfid); - fseek(restartfid, 0L, SEEK_SET); - - /*Allocate buffer: */ - femmodel_buffer=xNew(femmodel_size); - - /*Read buffer from file: */ - fread_return=fread(femmodel_buffer,femmodel_size,sizeof(char),restartfid); if(fread_return!=1)_error_("error reading the buffer from marshalled file!"); - femmodel_buffer_ini=femmodel_buffer; //keep track of the initial position, so as to free later. - - this->parameters->FindParam(&output_fid,OutputFilePointerEnum); - - /*Create new FemModel by demarshalling the buffer: */ - this->Marshall(&femmodel_buffer,NULL,MARSHALLING_BACKWARD); - - this->parameters->SetParam(output_fid,OutputFilePointerEnum); - - /*Reset position of buffer: */ - femmodel_buffer=femmodel_buffer_ini; - - /*Done, close file :*/ - pfclose(restartfid,restartfilename); - - /*Free ressources: */ - xDelete(restartfilename); - xDelete(femmodel_buffer); - -} -/*}}}*/ -void FemModel::Marshall(char** pmarshalled_data, int* pmarshalled_data_size, int marshall_direction){ /*{{{*/ - - int i; - int analysis_type; - - if(marshall_direction==MARSHALLING_BACKWARD){ - delete this->profiler; - delete this->loads; - delete this->materials; - delete this->parameters; - delete this->constraints; - delete this->results; - delete this->nodes; - delete this->vertices; - delete this->elements; - xDelete(this->analysis_type_list); - - this->profiler = new Profiler(); - this->loads = new Loads(); - this->materials = new Materials(); - this->parameters = new Parameters(); - this->constraints = new Constraints(); - this->results = new Results(); - this->nodes = new Nodes(); - this->vertices = new Vertices(); - this->elements = new Elements(); - } - - MARSHALLING_ENUM(FemModelEnum); - - MARSHALLING(solution_type); - MARSHALLING(analysis_counter); - MARSHALLING(nummodels); - MARSHALLING_DYNAMIC(analysis_type_list,int,nummodels); - - this->profiler->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); - this->loads->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); - this->materials->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); - this->parameters->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); - this->constraints->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); - this->results->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); - this->nodes->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); - this->vertices->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); - this->elements->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); - - if(marshall_direction==MARSHALLING_BACKWARD){ - //reset hooks for elements, loads and nodes: - this->elements->ResetHooks(); - this->loads->ResetHooks(); - this->materials->ResetHooks(); - - //do the post-processing of the datasets to get an FemModel that can actually run analyses: - for(i=0;ianalysis_type_list[i]; - SetCurrentConfiguration(analysis_type); - if(i==0) VerticesDofx(this->vertices,this->parameters); //only call once, we only have one set of vertices - SpcNodesx(this->nodes,this->constraints,this->parameters,analysis_type); - NodesDofx(this->nodes,this->parameters,analysis_type); - ConfigureObjectsx(this->elements,this->loads,this->nodes,this->vertices,this->materials,this->parameters); +void FemModel::Solve(void){/*{{{*/ + + /*profiling: */ + bool profiling = false; + IssmDouble solution_time; + IssmDouble solution_flops; + IssmDouble solution_memory; + + /*solution: */ + int solution_type; + void (*solutioncore)(FemModel*)=NULL; //core solution function pointer + + _printf0_("call computational core:\n"); + + /*Retrieve solution_type from parameters: */ + parameters->FindParam(&solution_type,SolutionTypeEnum); + + /*Figure out which solution core we are going to run with the current solution type: */ + WrapperCorePointerFromSolutionEnum(&solutioncore,this->parameters,solution_type); + + /*run solution core: */ + profiler->Tag(STARTCORE); + solutioncore(this); + profiler->Tag(FINISHCORE); + + /*run AD core if needed: */ + profiler->Tag(STARTADCORE); + ad_core(this); + profiler->Tag(FINISHADCORE); + + /*some profiling results for the core: */ + parameters->FindParam(&profiling,DebugProfilingEnum); + if(profiling){ + + solution_time=profiler->DeltaTime(STARTCORE,FINISHCORE); + solution_flops=profiler->DeltaFlops(STARTCORE,FINISHCORE); + solution_memory=profiler->Memory(FINISHCORE); + + _printf0_("Core solution elapsed time : " << solution_time << " Seconds\n"); + _printf0_("Core solution number of flops : " << solution_flops << " Flops\n"); + _printf0_("Core solution memory used : " << solution_memory << " Bytes\n"); + + /*Add to results: */ + results->AddObject(new GenericExternalResult(results->Size()+1, ProfilingSolutionTimeEnum, solution_time, 0, 0)); + results->AddObject(new GenericExternalResult(results->Size()+1, ProfilingCurrentMemEnum, solution_memory, 0, 0)); + results->AddObject(new GenericExternalResult(results->Size()+1, ProfilingCurrentFlopsEnum, solution_flops, 0, 0)); + + #ifdef _HAVE_ADOLC_ + solution_time=profiler->DeltaTime(STARTADCORE,FINISHADCORE); + solution_flops=profiler->DeltaFlops(STARTADCORE,FINISHADCORE); + solution_memory=profiler->Memory(FINISHADCORE); + + _printf0_("AD Solution elapsed time : " << solution_time << " Seconds\n"); + _printf0_("AD Solution number of flops : " << solution_flops << " Flops\n"); + _printf0_("AD Solution memory used : " << solution_memory << " Bytes\n"); + + #endif + + } +} +/*}}}*/ + +/*Modules:*/ +void FemModel::BalancethicknessMisfitx(IssmDouble* presponse){/*{{{*/ + + /*output: */ + IssmDouble J=0.; + IssmDouble J_sum; + + IssmDouble weight,vx,vy,H,dvx[2],dvy[2],dH[2]; + IssmDouble temp,Jdet,dhdt,groundedice_melting,surface_mass_balance; + IssmDouble* xyz_list = NULL; + IssmDouble dp[3]; + + /*Compute Misfit: */ + for(int i=0;iSize();i++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + + /*If on water, return 0: */ + if(!element->IsIceInElement()) continue; + + /* Get node coordinates*/ + element->GetVerticesCoordinates(&xyz_list); + Input* weights_input = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); + Input* surface_mass_balance_input = element->GetInput(SmbMassBalanceEnum); _assert_(surface_mass_balance_input); + Input* groundedice_melting_input = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(groundedice_melting_input); + Input* dhdt_input = element->GetInput(BalancethicknessThickeningRateEnum); _assert_(dhdt_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + /* Get Jacobian determinant: */ + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + + /*Get all parameters at gaussian point*/ + weights_input->GetInputValue(&weight,gauss,BalancethicknessMisfitEnum); + thickness_input->GetInputValue(&H, gauss); + thickness_input->GetInputDerivativeValue(&dH[0],xyz_list,gauss); + surface_mass_balance_input->GetInputValue(&surface_mass_balance,gauss); + groundedice_melting_input->GetInputValue(&groundedice_melting,gauss); + dhdt_input->GetInputValue(&dhdt,gauss); + vx_input->GetInputValue(&vx,gauss); + vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + vy_input->GetInputValue(&vy,gauss); + vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + + /*Balance thickness soft constraint J = 1/2 (div(Hv)-a)^2*/ + temp = vx*dH[0]+vy*dH[1]+H*(dvx[0]+dvy[1]) - (surface_mass_balance-groundedice_melting-dhdt); + J +=weight*1/2*temp*temp*Jdet*gauss->weight; } - //Reset current configuration: - analysis_type=this->analysis_type_list[analysis_counter]; - SetCurrentConfiguration(analysis_type); - } - -} -/*}}}*/ - -/*Modules:*/ -int FemModel::UpdateVertexPositionsx(void){ /*{{{*/ - - IssmDouble *surface = NULL; - IssmDouble *bed = NULL; - - if(VerboseSolution()) _printf0_(" updating vertices positions\n"); - - /*get vertex vectors for bed and thickness: */ - GetVectorFromInputsx(&surface ,this, SurfaceEnum,VertexPIdEnum); - GetVectorFromInputsx(&bed ,this, BaseEnum, VertexPIdEnum); - - /*Allocate vector*/ - Vector *vx=new Vector(vertices->NumberOfVertices()); - Vector *vy=new Vector(vertices->NumberOfVertices()); - Vector *vz=new Vector(vertices->NumberOfVertices()); - - /*Update verices new geometry: */ - for(int i=0;iSize();i++){ - Vertex* vertex=(Vertex*)vertices->GetObjectByOffset(i); - vertex->UpdatePosition(vx,vy,vz,parameters,surface,bed); - } - - /*Assemble mesh velocity*/ - vx->Assemble(); - vy->Assemble(); - vz->Assemble(); - - /*Update element inputs*/ - InputUpdateFromVectorx(this,vx,VxMeshEnum,VertexPIdEnum); - InputUpdateFromVectorx(this,vy,VyMeshEnum,VertexPIdEnum); - InputUpdateFromVectorx(this,vz,VzMeshEnum,VertexPIdEnum); - - /*Free ressources:*/ - delete vx; - delete vy; - delete vz; - xDelete(bed); - xDelete(surface); - return 1; -} -/*}}}*/ -void FemModel::UpdateConstraintsx(void){ /*{{{*/ - - IssmDouble time; - int analysis_type; - - /*retrieve parameters: */ - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - parameters->FindParam(&time,TimeEnum); - - /*start module: */ - if(VerboseModule()) _printf0_(" Updating constraints and active domain of analysis " << EnumToStringx(analysis_type) << " for time: " << time << "\n"); - - Analysis* analysis= EnumToAnalysis(analysis_type); - analysis->UpdateConstraints(this); - delete analysis; - - /*Second, constraints might be time dependent: */ - SpcNodesx(nodes,constraints,parameters,analysis_type); - - /*Now, update degrees of freedoms: */ - NodesDofx(nodes,parameters,analysis_type); - -} -/*}}}*/ + /*clean up and Return: */ + xDelete(xyz_list); + delete gauss; + } + + /*Sum all J from all cpus of the cluster:*/ + ISSM_MPI_Reduce (&J,&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + J=J_sum; + + /*Assign output pointers: */ + *presponse=J; + +}/*}}}*/ +void FemModel::CalvingRateDevx(){/*{{{*/ + + for(int i=0;iSize();i++){ + Element* element=dynamic_cast(this->elements->GetObjectByOffset(i)); + element->CalvingRateDev(); + } +} +/*}}}*/ +void FemModel::CalvingRateLevermannx(){/*{{{*/ + + for(int i=0;iSize();i++){ + Element* element=dynamic_cast(this->elements->GetObjectByOffset(i)); + element->CalvingRateLevermann(); + } +} +/*}}}*/ +void FemModel::CostFunctionx(IssmDouble* pJ,IssmDouble** pJlist,int* pn){/*{{{*/ + + /*Intermediary*/ + int num_responses; + int *responses = NULL; + Results *cost_functions = NULL; + + /*Recover parameters*/ + parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum); + parameters->FindParam(&responses,NULL,InversionCostFunctionsEnum); + + /*Get the value of all cost functions*/ + this->RequestedOutputsx(&cost_functions,responses,num_responses); + + /*Get and add all contributions one by one*/ + IssmDouble J=0; + IssmDouble* Jlist = xNew(num_responses); + for(int i=0;iGetObjectByOffset(i); + Jlist[i] = reCast(result->GetValue()); + J += Jlist[i]; + } + _assert_(cost_functions->Size()==num_responses); + + /*Assign output pointers: */ + delete cost_functions; + xDelete(responses); + if(pJ) *pJ = J; + if(pJlist) *pJlist = Jlist; + else xDelete(Jlist); + if(pn) *pn = num_responses; +} +/*}}}*/ +void FemModel::DeviatoricStressx(){/*{{{*/ + + for(int i=0;iSize();i++){ + Element* element=dynamic_cast(this->elements->GetObjectByOffset(i)); + element->ComputeDeviatoricStressTensor(); + } +} +/*}}}*/ +void FemModel::Divergencex(IssmDouble* pdiv){/*{{{*/ + + IssmDouble local_divergence=0; + IssmDouble total_divergence; + + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + local_divergence+=element->Divergence(); + } + ISSM_MPI_Reduce(&local_divergence,&total_divergence,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&total_divergence,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + + /*Assign output pointers: */ + *pdiv=total_divergence; + +}/*}}}*/ void FemModel::ElementOperationx(void (Element::*function)(void)){ /*{{{*/ @@ -869,4 +945,788 @@ } +} +/*}}}*/ +void FemModel::ElementResponsex(IssmDouble* presponse,int response_enum){/*{{{*/ + + int found=0; + int sumfound=0; + int cpu_found=-1; + int index; + IssmDouble response; + Element* element=NULL; + + /*retrieve element we are interested in: */ + this->parameters->FindParam(&index,IndexEnum); + int my_rank=IssmComm::GetRank(); + + /*now, go through our elements, and retrieve the one with this id: index: */ + for(int i=0;ielements->Size();i++){ + element=xDynamicCast(this->elements->GetObjectByOffset(i)); + if (element->Id()==index){ + found=1; + cpu_found=my_rank; + break; + } + } + + /*Broadcast whether we found the element: */ + ISSM_MPI_Allreduce ( &found,&sumfound,1,ISSM_MPI_INT,ISSM_MPI_SUM,IssmComm::GetComm()); + if(!sumfound)_error_("could not find material with id" << index << " to compute ElementResponse"); + + /*Ok, we found the element, compute responseocity: */ + if(my_rank==cpu_found){ + element->ElementResponse(&response,response_enum); + } + + /*Broadcast and plug into response: */ + ISSM_MPI_Allreduce ( &cpu_found,&cpu_found,1,ISSM_MPI_INT,ISSM_MPI_MAX,IssmComm::GetComm()); + ISSM_MPI_Bcast(&response,1,ISSM_MPI_DOUBLE,cpu_found,IssmComm::GetComm()); + + /*Assign output pointers: */ + *presponse=response; + +}/*}}}*/ +void FemModel::FloatingAreax(IssmDouble* pV){/*{{{*/ + + IssmDouble local_floating_area= 0; + IssmDouble total_floating_area; + + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + local_floating_area+=element->FloatingArea(); + } + ISSM_MPI_Reduce(&local_floating_area,&total_floating_area,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&total_floating_area,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + + /*Assign output pointers: */ + *pV=total_floating_area; + +}/*}}}*/ +void FemModel::GetInputLocalMinMaxOnNodesx(IssmDouble** pmin,IssmDouble** pmax,IssmDouble* ug){/*{{{*/ + + /*Get vector sizes for current configuration*/ + int configuration_type; + this->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + int numnodes = this->nodes->NumberOfNodes(configuration_type); + + /*Initialize output vectors*/ + IssmDouble* uLmin_local = xNew(numnodes); + IssmDouble* uLmax_local = xNew(numnodes); + IssmDouble* uLmin = xNew(numnodes); + IssmDouble* uLmax = xNew(numnodes); + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + element->GetInputLocalMinMaxOnNodes(uLmin_local,uLmax_local,ug); + } + + /*Synchronize all CPUs*/ + ISSM_MPI_Allreduce((void*)uLmin_local,(void*)uLmin,numnodes,ISSM_MPI_DOUBLE,ISSM_MPI_MIN,IssmComm::GetComm()); + ISSM_MPI_Allreduce((void*)uLmax_local,(void*)uLmax,numnodes,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,IssmComm::GetComm()); + xDelete(uLmin_local); + xDelete(uLmax_local); + + /*Assign output pointers: */ + *pmin=uLmin; + *pmax=uLmax; + +}/*}}}*/ +void FemModel::GroundedAreax(IssmDouble* pV){/*{{{*/ + + IssmDouble local_grounded_area= 0; + IssmDouble total_grounded_area; + + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + local_grounded_area+=element->GroundedArea(); + } + ISSM_MPI_Reduce(&local_grounded_area,&total_grounded_area,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&total_grounded_area,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + + /*Assign output pointers: */ + *pV=total_grounded_area; + +}/*}}}*/ +void FemModel::IceMassx(IssmDouble* pM){/*{{{*/ + + IssmDouble local_ice_mass = 0; + IssmDouble total_ice_mass; + + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + local_ice_mass+=element->IceMass(); + } + ISSM_MPI_Reduce(&local_ice_mass,&total_ice_mass,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&total_ice_mass,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + + /*Assign output pointers: */ + *pM=total_ice_mass; + +}/*}}}*/ +void FemModel::IceVolumeAboveFloatationx(IssmDouble* pV){/*{{{*/ + + IssmDouble local_ice_volume_af = 0; + IssmDouble total_ice_volume_af; + + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + local_ice_volume_af+=element->IceVolumeAboveFloatation(); + } + ISSM_MPI_Reduce(&local_ice_volume_af,&total_ice_volume_af,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&total_ice_volume_af,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + + /*Assign output pointers: */ + *pV=total_ice_volume_af; + +}/*}}}*/ +void FemModel::IceVolumex(IssmDouble* pV){/*{{{*/ + + IssmDouble local_ice_volume = 0; + IssmDouble total_ice_volume; + + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + local_ice_volume+=element->IceVolume(); + } + ISSM_MPI_Reduce(&local_ice_volume,&total_ice_volume,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&total_ice_volume,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + + /*Assign output pointers: */ + *pV=total_ice_volume; + +}/*}}}*/ +void FemModel::MassFluxx(IssmDouble* pmass_flux){/*{{{*/ + + int i,j; + Element *element = NULL; + int element_id; + bool ispresent = false; + IssmDouble mass_flux = 0; + IssmDouble all_mass_flux = 0; + int counter; + IssmDouble **array = NULL; + int M; + int *mdims_array = NULL; + int *ndims_array = NULL; + IssmDouble *segments = NULL; + int num_segments; + + /*First, figure out which segment to compute our mass flux on. Start with retrieving qmu_mass_flux_segments: */ + this->parameters->FindParam(&ispresent,MassFluxSegmentsPresentEnum); + if(!ispresent)_error_("no mass flux segments available!"); + this->parameters->FindParam(&array,&M,&mdims_array,&ndims_array,MassFluxSegmentsEnum); + + /*Retrieve index of segments being used for MassFlux computation: */ + parameters->FindParam(&counter,IndexEnum); + + /*retrieve segments from array: */ + segments = array[counter-1]; //matlab to "C" indexing + num_segments = mdims_array[counter-1]; + + /*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;i(*(segments+5*i+4)); + for(j=0;jSize();j++){ + element=(Element*)this->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); + break; + } + } + } + + ISSM_MPI_Allreduce ( (void*)&mass_flux,(void*)&all_mass_flux,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); + mass_flux=all_mass_flux; + + /*Free ressources:*/ + for(i=0;i(matrix); + } + xDelete(mdims_array); + xDelete(ndims_array); + xDelete(array); + + /*Assign output pointers: */ + *pmass_flux=mass_flux; + +}/*}}}*/ +void FemModel::MaxAbsVxx(IssmDouble* pmaxabsvx){/*{{{*/ + + int i; + IssmDouble maxabsvx; + IssmDouble node_maxabsvx; + IssmDouble element_maxabsvx; + + /*Go through elements, and request velocity: */ + maxabsvx=-INFINITY; + for(i=0;ielements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + element_maxabsvx=element->inputs->MaxAbs(VxEnum); + if(element_maxabsvx>maxabsvx) maxabsvx=element_maxabsvx; + } + + /*Figure out maximum across the cluster: */ + ISSM_MPI_Reduce(&maxabsvx,&node_maxabsvx,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&node_maxabsvx,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + maxabsvx=node_maxabsvx; + + /*Assign output pointers:*/ + *pmaxabsvx=maxabsvx; + +}/*}}}*/ +void FemModel::MaxAbsVyx(IssmDouble* pmaxabsvy){/*{{{*/ + + int i; + IssmDouble maxabsvy; + IssmDouble node_maxabsvy; + IssmDouble element_maxabsvy; + + /*Go through elements, and request velocity: */ + maxabsvy=-INFINITY; + for(i=0;ielements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + element_maxabsvy=element->inputs->MaxAbs(VyEnum); + if(element_maxabsvy>maxabsvy) maxabsvy=element_maxabsvy; + } + + /*Figure out maximum across the cluster: */ + ISSM_MPI_Reduce(&maxabsvy,&node_maxabsvy,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&node_maxabsvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + maxabsvy=node_maxabsvy; + + /*Assign output pointers:*/ + *pmaxabsvy=maxabsvy; + +}/*}}}*/ +void FemModel::MaxAbsVzx(IssmDouble* pmaxabsvz){/*{{{*/ + + int i; + IssmDouble maxabsvz; + IssmDouble node_maxabsvz; + IssmDouble element_maxabsvz; + + /*Go through elements, and request velocity: */ + maxabsvz=-INFINITY; + for(i=0;ielements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + element_maxabsvz=element->inputs->MaxAbs(VzEnum); + if(element_maxabsvz>maxabsvz) maxabsvz=element_maxabsvz; + } + + /*Figure out maximum across the cluster: */ + ISSM_MPI_Reduce(&maxabsvz,&node_maxabsvz,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&node_maxabsvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + maxabsvz=node_maxabsvz; + + /*Assign output pointers:*/ + *pmaxabsvz=maxabsvz; + +}/*}}}*/ +void FemModel::MaxDivergencex(IssmDouble* pdiv){/*{{{*/ + + IssmDouble local_divergence; + IssmDouble node_max_divergence; + IssmDouble max_divergence = -INFINITY; + + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + local_divergence=element->Divergence(); + if(fabs(local_divergence)>max_divergence) max_divergence=fabs(local_divergence); + } + ISSM_MPI_Reduce(&max_divergence,&node_max_divergence,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&node_max_divergence,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + max_divergence=node_max_divergence; + + /*Assign output pointers: */ + *pdiv=max_divergence; + +}/*}}}*/ +void FemModel::MaxVelx(IssmDouble* pmaxvel){/*{{{*/ + + int i; + IssmDouble maxvel; + IssmDouble node_maxvel; + IssmDouble element_maxvel; + + /*Go through elements, and request velocity: */ + maxvel=-INFINITY; + for(i=0;ielements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + element_maxvel = element->inputs->Max(VelEnum); + if(element_maxvel>maxvel) maxvel=element_maxvel; + } + + /*Figure out maximum across the cluster: */ + ISSM_MPI_Reduce(&maxvel,&node_maxvel,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&node_maxvel,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + maxvel=node_maxvel; + + /*Assign output pointers:*/ + *pmaxvel=maxvel; + +}/*}}}*/ +void FemModel::MaxVxx(IssmDouble* pmaxvx){/*{{{*/ + + int i; + IssmDouble maxvx; + IssmDouble node_maxvx; + IssmDouble element_maxvx; + + /*Go through elements, and request velocity: */ + maxvx=-INFINITY; + for(i=0;ielements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + element_maxvx = element->inputs->Max(VxEnum); + if(element_maxvx>maxvx) maxvx=element_maxvx; + } + + /*Figure out maximum across the cluster: */ + ISSM_MPI_Reduce(&maxvx,&node_maxvx,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&node_maxvx,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + maxvx=node_maxvx; + + /*Assign output pointers:*/ + *pmaxvx=maxvx; + +}/*}}}*/ +void FemModel::MaxVyx(IssmDouble* pmaxvy){/*{{{*/ + + int i; + IssmDouble maxvy; + IssmDouble node_maxvy; + IssmDouble element_maxvy; + + /*Go through elements, and request velocity: */ + maxvy=-INFINITY; + for(i=0;ielements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + element_maxvy = element->inputs->Max(VyEnum); + if(element_maxvy>maxvy) maxvy=element_maxvy; + } + + /*Figure out maximum across the cluster: */ + ISSM_MPI_Reduce(&maxvy,&node_maxvy,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&node_maxvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + maxvy=node_maxvy; + + /*Assign output pointers:*/ + *pmaxvy=maxvy; + +}/*}}}*/ +void FemModel::MaxVzx(IssmDouble* pmaxvz){/*{{{*/ + + int i; + IssmDouble maxvz; + IssmDouble node_maxvz; + IssmDouble element_maxvz; + + /*Go through elements, and request velocity: */ + maxvz=-INFINITY; + for(i=0;ielements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + element_maxvz = element->inputs->Max(VzEnum); + if(element_maxvz>maxvz) maxvz=element_maxvz; + } + + /*Figure out maximum across the cluster: */ + ISSM_MPI_Reduce(&maxvz,&node_maxvz,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&node_maxvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + maxvz=node_maxvz; + + /*Assign output pointers:*/ + *pmaxvz=maxvz; + +}/*}}}*/ +void FemModel::MinVelx(IssmDouble* pminvel){/*{{{*/ + + int i; + IssmDouble minvel; + IssmDouble node_minvel; + IssmDouble element_minvel; + + /*Go through elements, and request velocity: */ + minvel=INFINITY; + for(i=0;ielements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + element_minvel = element->inputs->Min(VelEnum); + if(element_minvelelements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + element_minvx = element->inputs->Min(VxEnum); + if(element_minvxelements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + element_minvy = element->inputs->Min(VyEnum); + if(element_minvyelements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + element_minvz = element->inputs->Min(VzEnum); + if(element_minvzparameters->FindParam(&num_controls,InversionNumControlParametersEnum); + this->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum); + this->parameters->FindParam(&step,StepEnum); + this->parameters->FindParam(&time,TimeEnum); + + for(int i=0;i *vector_control = new Vector(this->vertices->NumberOfVertices()); + Vector *vector_gradient = new Vector(this->vertices->NumberOfVertices()); + + /*Fill in vector*/ + for(int j=0;jSize();j++){ + Element* element=(Element*)elements->GetObjectByOffset(j); + element->ControlToVectors(vector_control,vector_gradient,control_enum); + } + vector_control->Assemble(); + vector_gradient->Assemble(); + + results->AddResult(new GenericExternalResult*>(results->Size()+1,control_enum,vector_control ,step,time)); + results->AddResult(new GenericExternalResult*>(results->Size()+1,gradient_enum,vector_gradient,step,time)); + } + + /*Clean up and return*/ + xDelete(control_type); +} +/*}}}*/ +void FemModel::RequestedDependentsx(void){/*{{{*/ + + bool isautodiff = false; + IssmDouble output_value; + + int num_dependents; + IssmPDouble *dependents; + DataSet* dependent_objects=NULL; + int my_rank=IssmComm::GetRank(); + + /*AD mode on?: */ + parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum); + + if(isautodiff){ + #ifdef _HAVE_ADOLC_ + parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum); + parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum); + if(num_dependents){ + dependents=xNew(num_dependents); + + /*Go through our dependent variables, and compute the response:*/ + for(int i=0;iSize();i++){ + DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i); + dep->Responsex(&output_value,this); + if (my_rank==0) { + output_value>>=dependents[i]; + } + } + } + delete dependent_objects; + if(num_dependents)xDelete(dependents); + #else + _error_("Should not be requesting dependents when an AD library is not available!"); + #endif + } +} +/*}}}*/ +void FemModel::RequestedOutputsx(Results **presults,char** requested_outputs, int numoutputs, bool save_results){/*{{{*/ + + /*Intermediaries*/ + bool isvec,results_on_nodes; + int step,output_enum; + IssmDouble time; + IssmDouble double_result; + const char *output_string = NULL; + + /*recover results*/ + Results* results = *presults; + if(!results) results = new Results(); + + /*Get time and step*/ + parameters->FindParam(&step,StepEnum); + parameters->FindParam(&time,TimeEnum); + parameters->FindParam(&results_on_nodes,SettingsResultsOnNodesEnum); + + /*Go through all requested output*/ + for(int i=0;iAddResult(new GenericExternalResult(results->Size()+1,output_string,reCast(double_result),step,time)); + continue; + } + } + else{ + /*last chance for the output definition, if the enum is one of Outputdefinition[1-10]Enum:*/ + if(output_enum>=Outputdefinition1Enum && output_enum <=Outputdefinition10Enum){ + double_result = OutputDefinitionsResponsex(this,output_enum); + if(save_results){ + results->AddResult(new GenericExternalResult(results->Size()+1,output_string,reCast(double_result),step,time)); + continue; + } + } + else{ + switch(output_enum){ + + /*Scalar output*/ + case DivergenceEnum: this->Divergencex(&double_result); break; + case MaxDivergenceEnum: this->MaxDivergencex(&double_result); break; + case IceMassEnum: this->IceMassx(&double_result); break; + case IceVolumeEnum: this->IceVolumex(&double_result); break; + case IceVolumeAboveFloatationEnum: this->IceVolumeAboveFloatationx(&double_result); break; + case GroundedAreaEnum: this->GroundedAreax(&double_result); break; + case FloatingAreaEnum: this->FloatingAreax(&double_result); break; + case MinVelEnum: this->MinVelx(&double_result); break; + case MaxVelEnum: this->MaxVelx(&double_result); break; + case MinVxEnum: this->MinVxx(&double_result); break; + case MaxVxEnum: this->MaxVxx(&double_result); break; + case MaxAbsVxEnum: this->MaxAbsVxx(&double_result); break; + case MinVyEnum: this->MinVyx(&double_result); break; + case MaxVyEnum: this->MaxVyx(&double_result); break; + case MaxAbsVyEnum: this->MaxAbsVyx(&double_result); break; + case MinVzEnum: this->MinVzx(&double_result); break; + case MaxVzEnum: this->MaxVzx(&double_result); break; + case MaxAbsVzEnum: this->MaxAbsVzx(&double_result); break; + case MassFluxEnum: this->MassFluxx(&double_result); break; + case TotalFloatingBmbEnum: this->TotalFloatingBmbx(&double_result); break; + case TotalGroundedBmbEnum: this->TotalGroundedBmbx(&double_result); break; + case TotalSmbEnum: this->TotalSmbx(&double_result); break; + + /*Scalar control output*/ + case SurfaceAbsVelMisfitEnum: SurfaceAbsVelMisfitx(&double_result,elements,nodes,vertices,loads,materials,parameters); break; + case SurfaceRelVelMisfitEnum: SurfaceRelVelMisfitx(&double_result,elements,nodes,vertices,loads,materials,parameters); break; + case SurfaceLogVelMisfitEnum: SurfaceLogVelMisfitx(&double_result,elements,nodes,vertices,loads,materials,parameters); break; + case SurfaceLogVxVyMisfitEnum: SurfaceLogVxVyMisfitx(&double_result,elements,nodes,vertices,loads,materials,parameters); break; + case SurfaceAverageVelMisfitEnum: SurfaceAverageVelMisfitx(&double_result,this); break; + case ThicknessAbsMisfitEnum: ThicknessAbsMisfitx(&double_result,elements,nodes,vertices,loads,materials,parameters); break; + case ThicknessAbsGradientEnum: this->ThicknessAbsGradientx(&double_result); break; + case ThicknessAlongGradientEnum: ThicknessAlongGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters); break; + case ThicknessAcrossGradientEnum: ThicknessAcrossGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters); break; + case ThicknessPositiveEnum: this->ThicknessPositivex(&double_result); break; + case RheologyBbarAbsGradientEnum: RheologyBbarAbsGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters); break; + case RheologyBAbsGradientEnum: RheologyBAbsGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters); break; + case DragCoefficientAbsGradientEnum:DragCoefficientAbsGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters); break; + case BalancethicknessMisfitEnum: BalancethicknessMisfitx(&double_result); break; + case SurfaceAbsMisfitEnum: SurfaceAbsMisfitx(&double_result); break; + + /*Vector */ + default: + + /*Vector layout*/ + int interpolation,nodesperelement,size,nlines,ncols,array_size; + int rank_interpolation=-1,rank_nodesperelement=-1,rank_arraysize=-1,max_rank_arraysize=0; + bool isarray=false; + + /*Get interpolation (and compute input if necessary)*/ + for(int j=0;jSize();j++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(j)); + element->ResultInterpolation(&rank_interpolation,&rank_nodesperelement,&rank_arraysize,output_enum); + if (rank_arraysize>max_rank_arraysize)max_rank_arraysize=rank_arraysize; + } + rank_arraysize=max_rank_arraysize; + + /*Broadcast for cpus that do not have any elements*/ + ISSM_MPI_Reduce(&rank_interpolation,&interpolation,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm()); + ISSM_MPI_Reduce(&rank_nodesperelement,&nodesperelement,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm()); + ISSM_MPI_Reduce(&rank_arraysize,&array_size,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&interpolation,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&nodesperelement,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&array_size,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + + if(results_on_nodes){ + + /*Allocate matrices*/ + int nbe = this->elements->NumberOfElements(); + IssmDouble* values = xNewZeroInit(nbe*nodesperelement); + IssmDouble* allvalues = xNew(nbe*nodesperelement); + + /*Fill-in matrix*/ + for(int j=0;jSize();j++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(j)); + element->ResultToPatch(values,nodesperelement,output_enum); + } + + /*Gather from all cpus*/ + ISSM_MPI_Allreduce((void*)values,(void*)allvalues,nbe*nodesperelement,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); + xDelete(values); + + if(save_results)results->AddResult(new GenericExternalResult(results->Size()+1,output_enum,allvalues,nbe,nodesperelement,step,time)); + xDelete(allvalues); + + } + else{ + + /*Allocate vector depending on interpolation*/ + switch(interpolation){ + case P0Enum: isarray = false; size = this->elements->NumberOfElements(); break; + case P1Enum: isarray = false; size = this->vertices->NumberOfVertices(); break; + case P0ArrayEnum: isarray = true; nlines = this->elements->NumberOfElements(); ncols= array_size; break; + default: _error_("Interpolation "< *vector_result = new Vector(size); + + /*Fill in vector*/ + for(int j=0;jSize();j++){ + Element* element=(Element*)elements->GetObjectByOffset(j); + element->ResultToVector(vector_result,output_enum); + } + vector_result->Assemble(); + + if(save_results)results->AddResult(new GenericExternalResult*>(results->Size()+1,output_enum,vector_result,step,time)); + } + else{ + IssmDouble* values = xNewZeroInit(nlines*ncols); + IssmDouble* allvalues = xNew(nlines*ncols); + + /*Fill-in matrix*/ + for(int j=0;jSize();j++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(j)); + element->ResultToMatrix(values,ncols, output_enum); + } + /*Gather from all cpus*/ + ISSM_MPI_Allreduce((void*)values,(void*)allvalues,ncols*nlines,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); + xDelete(values); + + if(save_results)results->AddResult(new GenericExternalResult(results->Size()+1,output_enum,allvalues,nlines,ncols,step,time)); + xDelete(allvalues); + } + } + isvec = true; + break; + } + } + + } + + /*Add result to Results*/ + if(!isvec && save_results){ + results->AddResult(new GenericExternalResult(results->Size()+1,output_string,reCast(double_result),step,time)); + } + } + + /*Assign pointer and clean up*/ + *presults = results; +} +/*}}}*/ +void FemModel::RequestedOutputsx(Results **presults,int* requested_outputs, int numoutputs,bool save_results){/*{{{*/ + + /*Convert list of enums to list of string*/ + char** enumlist = xNew(numoutputs); + for(int i=0;iRequestedOutputsx(presults,enumlist,numoutputs,save_results); + + /*clean up and return*/ + for(int i=0;i(enumlist[i]); + xDelete(enumlist); + return; } /*}}}*/ @@ -975,919 +1835,29 @@ } /*}}}*/ -void FemModel::RequestedOutputsx(Results **presults,char** requested_outputs, int numoutputs, bool save_results){/*{{{*/ - - /*Intermediaries*/ - bool isvec,results_on_nodes; - int step,output_enum; - IssmDouble time; - IssmDouble double_result; - const char *output_string = NULL; - - /*recover results*/ - Results* results = *presults; - if(!results) results = new Results(); - - /*Get time and step*/ - parameters->FindParam(&step,StepEnum); - parameters->FindParam(&time,TimeEnum); - parameters->FindParam(&results_on_nodes,SettingsResultsOnNodesEnum); - - /*Go through all requested output*/ - for(int i=0;iAddResult(new GenericExternalResult(results->Size()+1,output_string,reCast(double_result),step,time)); - continue; - } - } - else{ - /*last chance for the output definition, if the enum is one of Outputdefinition[1-10]Enum:*/ - if(output_enum>=Outputdefinition1Enum && output_enum <=Outputdefinition10Enum){ - double_result = OutputDefinitionsResponsex(this,output_enum); - if(save_results){ - results->AddResult(new GenericExternalResult(results->Size()+1,output_string,reCast(double_result),step,time)); - continue; - } - } - else{ - switch(output_enum){ - - /*Scalar output*/ - case DivergenceEnum: this->Divergencex(&double_result); break; - case MaxDivergenceEnum: this->MaxDivergencex(&double_result); break; - case IceMassEnum: this->IceMassx(&double_result); break; - case IceVolumeEnum: this->IceVolumex(&double_result); break; - case IceVolumeAboveFloatationEnum: this->IceVolumeAboveFloatationx(&double_result); break; - case GroundedAreaEnum: this->GroundedAreax(&double_result); break; - case FloatingAreaEnum: this->FloatingAreax(&double_result); break; - case MinVelEnum: this->MinVelx(&double_result); break; - case MaxVelEnum: this->MaxVelx(&double_result); break; - case MinVxEnum: this->MinVxx(&double_result); break; - case MaxVxEnum: this->MaxVxx(&double_result); break; - case MaxAbsVxEnum: this->MaxAbsVxx(&double_result); break; - case MinVyEnum: this->MinVyx(&double_result); break; - case MaxVyEnum: this->MaxVyx(&double_result); break; - case MaxAbsVyEnum: this->MaxAbsVyx(&double_result); break; - case MinVzEnum: this->MinVzx(&double_result); break; - case MaxVzEnum: this->MaxVzx(&double_result); break; - case MaxAbsVzEnum: this->MaxAbsVzx(&double_result); break; - case MassFluxEnum: this->MassFluxx(&double_result); break; - case TotalFloatingBmbEnum: this->TotalFloatingBmbx(&double_result); break; - case TotalGroundedBmbEnum: this->TotalGroundedBmbx(&double_result); break; - case TotalSmbEnum: this->TotalSmbx(&double_result); break; - - /*Scalar control output*/ - case SurfaceAbsVelMisfitEnum: SurfaceAbsVelMisfitx(&double_result,elements,nodes,vertices,loads,materials,parameters); break; - case SurfaceRelVelMisfitEnum: SurfaceRelVelMisfitx(&double_result,elements,nodes,vertices,loads,materials,parameters); break; - case SurfaceLogVelMisfitEnum: SurfaceLogVelMisfitx(&double_result,elements,nodes,vertices,loads,materials,parameters); break; - case SurfaceLogVxVyMisfitEnum: SurfaceLogVxVyMisfitx(&double_result,elements,nodes,vertices,loads,materials,parameters); break; - case SurfaceAverageVelMisfitEnum: SurfaceAverageVelMisfitx(&double_result,this); break; - case ThicknessAbsMisfitEnum: ThicknessAbsMisfitx(&double_result,elements,nodes,vertices,loads,materials,parameters); break; - case ThicknessAbsGradientEnum: this->ThicknessAbsGradientx(&double_result); break; - case ThicknessAlongGradientEnum: ThicknessAlongGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters); break; - case ThicknessAcrossGradientEnum: ThicknessAcrossGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters); break; - case ThicknessPositiveEnum: this->ThicknessPositivex(&double_result); break; - case RheologyBbarAbsGradientEnum: RheologyBbarAbsGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters); break; - case RheologyBAbsGradientEnum: RheologyBAbsGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters); break; - case DragCoefficientAbsGradientEnum:DragCoefficientAbsGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters); break; - case BalancethicknessMisfitEnum: BalancethicknessMisfitx(&double_result); break; - case SurfaceAbsMisfitEnum: SurfaceAbsMisfitx(&double_result); break; - - /*Vector */ - default: - - /*Vector layout*/ - int interpolation,nodesperelement,size,nlines,ncols,array_size; - int rank_interpolation=-1,rank_nodesperelement=-1,rank_arraysize=-1,max_rank_arraysize=0; - bool isarray=false; - - /*Get interpolation (and compute input if necessary)*/ - for(int j=0;jSize();j++){ - Element* element=xDynamicCast(this->elements->GetObjectByOffset(j)); - element->ResultInterpolation(&rank_interpolation,&rank_nodesperelement,&rank_arraysize,output_enum); - if (rank_arraysize>max_rank_arraysize)max_rank_arraysize=rank_arraysize; - } - rank_arraysize=max_rank_arraysize; - - /*Broadcast for cpus that do not have any elements*/ - ISSM_MPI_Reduce(&rank_interpolation,&interpolation,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm()); - ISSM_MPI_Reduce(&rank_nodesperelement,&nodesperelement,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm()); - ISSM_MPI_Reduce(&rank_arraysize,&array_size,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(&interpolation,1,ISSM_MPI_INT,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(&nodesperelement,1,ISSM_MPI_INT,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(&array_size,1,ISSM_MPI_INT,0,IssmComm::GetComm()); - - if(results_on_nodes){ - - /*Allocate matrices*/ - int nbe = this->elements->NumberOfElements(); - IssmDouble* values = xNewZeroInit(nbe*nodesperelement); - IssmDouble* allvalues = xNew(nbe*nodesperelement); - - /*Fill-in matrix*/ - for(int j=0;jSize();j++){ - Element* element=xDynamicCast(this->elements->GetObjectByOffset(j)); - element->ResultToPatch(values,nodesperelement,output_enum); - } - - /*Gather from all cpus*/ - ISSM_MPI_Allreduce((void*)values,(void*)allvalues,nbe*nodesperelement,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); - xDelete(values); - - if(save_results)results->AddResult(new GenericExternalResult(results->Size()+1,output_enum,allvalues,nbe,nodesperelement,step,time)); - xDelete(allvalues); - - } - else{ - - /*Allocate vector depending on interpolation*/ - switch(interpolation){ - case P0Enum: isarray = false; size = this->elements->NumberOfElements(); break; - case P1Enum: isarray = false; size = this->vertices->NumberOfVertices(); break; - case P0ArrayEnum: isarray = true; nlines = this->elements->NumberOfElements(); ncols= array_size; break; - default: _error_("Interpolation "< *vector_result = new Vector(size); - - /*Fill in vector*/ - for(int j=0;jSize();j++){ - Element* element=(Element*)elements->GetObjectByOffset(j); - element->ResultToVector(vector_result,output_enum); - } - vector_result->Assemble(); - - if(save_results)results->AddResult(new GenericExternalResult*>(results->Size()+1,output_enum,vector_result,step,time)); - } - else{ - IssmDouble* values = xNewZeroInit(nlines*ncols); - IssmDouble* allvalues = xNew(nlines*ncols); - - /*Fill-in matrix*/ - for(int j=0;jSize();j++){ - Element* element=xDynamicCast(this->elements->GetObjectByOffset(j)); - element->ResultToMatrix(values,ncols, output_enum); - } - /*Gather from all cpus*/ - ISSM_MPI_Allreduce((void*)values,(void*)allvalues,ncols*nlines,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); - xDelete(values); - - if(save_results)results->AddResult(new GenericExternalResult(results->Size()+1,output_enum,allvalues,nlines,ncols,step,time)); - xDelete(allvalues); - } - } - isvec = true; - break; - } - } - - } - - /*Add result to Results*/ - if(!isvec && save_results){ - results->AddResult(new GenericExternalResult(results->Size()+1,output_string,reCast(double_result),step,time)); - } - } - - /*Assign pointer and clean up*/ - *presults = results; -} -/*}}}*/ -void FemModel::RequestedOutputsx(Results **presults,int* requested_outputs, int numoutputs,bool save_results){/*{{{*/ - - /*Convert list of enums to list of string*/ - char** enumlist = xNew(numoutputs); - for(int i=0;iRequestedOutputsx(presults,enumlist,numoutputs,save_results); - - /*clean up and return*/ - for(int i=0;i(enumlist[i]); - xDelete(enumlist); - return; -} -/*}}}*/ -void FemModel::RequestedDependentsx(void){/*{{{*/ - - bool isautodiff = false; - IssmDouble output_value; - - int num_dependents; - IssmPDouble *dependents; - DataSet* dependent_objects=NULL; - int my_rank=IssmComm::GetRank(); - - /*AD mode on?: */ - parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum); - - if(isautodiff){ - #ifdef _HAVE_ADOLC_ - parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum); - parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum); - if(num_dependents){ - dependents=xNew(num_dependents); - - /*Go through our dependent variables, and compute the response:*/ - for(int i=0;iSize();i++){ - DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i); - dep->Responsex(&output_value,this); - if (my_rank==0) { - output_value>>=dependents[i]; - } - } - } - delete dependent_objects; - if(num_dependents)xDelete(dependents); - #else - _error_("Should not be requesting dependents when an AD library is not available!"); - #endif - } -} -/*}}}*/ -void FemModel::TimeAdaptx(IssmDouble* pdt){/*{{{*/ - - int i; - - /*output: */ - IssmDouble dt; - - /*intermediary: */ - Element *element = NULL; - IssmDouble min_dt = 0; - IssmDouble node_min_dt = 0; - - /*Go through elements, and figure out the minimum of the time steps for each element (using CFL criterion): */ - element=(Element*)elements->GetObjectByOffset(0); min_dt=element->TimeAdapt(); - - for (i=1;iSize();i++){ - element=xDynamicCast(elements->GetObjectByOffset(i)); - dt=element->TimeAdapt(); - if(dtSize();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->UpdateConstraintsExtrudeFromBase(); - } - -} -/*}}}*/ -void FemModel::UpdateConstraintsExtrudeFromTopx(void){ /*{{{*/ + Element* element=dynamic_cast(this->elements->GetObjectByOffset(i)); + element->StrainRateparallel(); + } +} +/*}}}*/ +void FemModel::StrainRateperpendicularx(){/*{{{*/ for(int i=0;iSize();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->UpdateConstraintsExtrudeFromTop(); - } - -} -/*}}}*/ -void FemModel::MassFluxx(IssmDouble* pmass_flux){/*{{{*/ - - int i,j; - Element *element = NULL; - int element_id; - bool ispresent = false; - IssmDouble mass_flux = 0; - IssmDouble all_mass_flux = 0; - int counter; - IssmDouble **array = NULL; - int M; - int *mdims_array = NULL; - int *ndims_array = NULL; - IssmDouble *segments = NULL; - int num_segments; - - /*First, figure out which segment to compute our mass flux on. Start with retrieving qmu_mass_flux_segments: */ - this->parameters->FindParam(&ispresent,MassFluxSegmentsPresentEnum); - if(!ispresent)_error_("no mass flux segments available!"); - this->parameters->FindParam(&array,&M,&mdims_array,&ndims_array,MassFluxSegmentsEnum); - - /*Retrieve index of segments being used for MassFlux computation: */ - parameters->FindParam(&counter,IndexEnum); - - /*retrieve segments from array: */ - segments = array[counter-1]; //matlab to "C" indexing - num_segments = mdims_array[counter-1]; - - /*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;i(*(segments+5*i+4)); - for(j=0;jSize();j++){ - element=(Element*)this->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); - break; - } - } - } - - ISSM_MPI_Allreduce ( (void*)&mass_flux,(void*)&all_mass_flux,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); - mass_flux=all_mass_flux; - - /*Free ressources:*/ - for(i=0;i(matrix); - } - xDelete(mdims_array); - xDelete(ndims_array); - xDelete(array); - - /*Assign output pointers: */ - *pmass_flux=mass_flux; - -}/*}}}*/ -void FemModel::MaxAbsVxx(IssmDouble* pmaxabsvx){/*{{{*/ - - int i; - IssmDouble maxabsvx; - IssmDouble node_maxabsvx; - IssmDouble element_maxabsvx; - - /*Go through elements, and request velocity: */ - maxabsvx=-INFINITY; - for(i=0;ielements->Size();i++){ - Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - element_maxabsvx=element->inputs->MaxAbs(VxEnum); - if(element_maxabsvx>maxabsvx) maxabsvx=element_maxabsvx; - } - - /*Figure out maximum across the cluster: */ - ISSM_MPI_Reduce(&maxabsvx,&node_maxabsvx,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&node_maxabsvx,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - maxabsvx=node_maxabsvx; - - /*Assign output pointers:*/ - *pmaxabsvx=maxabsvx; - -}/*}}}*/ -void FemModel::MaxAbsVyx(IssmDouble* pmaxabsvy){/*{{{*/ - - int i; - IssmDouble maxabsvy; - IssmDouble node_maxabsvy; - IssmDouble element_maxabsvy; - - /*Go through elements, and request velocity: */ - maxabsvy=-INFINITY; - for(i=0;ielements->Size();i++){ - Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - element_maxabsvy=element->inputs->MaxAbs(VyEnum); - if(element_maxabsvy>maxabsvy) maxabsvy=element_maxabsvy; - } - - /*Figure out maximum across the cluster: */ - ISSM_MPI_Reduce(&maxabsvy,&node_maxabsvy,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&node_maxabsvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - maxabsvy=node_maxabsvy; - - /*Assign output pointers:*/ - *pmaxabsvy=maxabsvy; - -}/*}}}*/ -void FemModel::MaxAbsVzx(IssmDouble* pmaxabsvz){/*{{{*/ - - int i; - IssmDouble maxabsvz; - IssmDouble node_maxabsvz; - IssmDouble element_maxabsvz; - - /*Go through elements, and request velocity: */ - maxabsvz=-INFINITY; - for(i=0;ielements->Size();i++){ - Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - element_maxabsvz=element->inputs->MaxAbs(VzEnum); - if(element_maxabsvz>maxabsvz) maxabsvz=element_maxabsvz; - } - - /*Figure out maximum across the cluster: */ - ISSM_MPI_Reduce(&maxabsvz,&node_maxabsvz,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&node_maxabsvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - maxabsvz=node_maxabsvz; - - /*Assign output pointers:*/ - *pmaxabsvz=maxabsvz; - -}/*}}}*/ -void FemModel::MaxVelx(IssmDouble* pmaxvel){/*{{{*/ - - int i; - IssmDouble maxvel; - IssmDouble node_maxvel; - IssmDouble element_maxvel; - - /*Go through elements, and request velocity: */ - maxvel=-INFINITY; - for(i=0;ielements->Size();i++){ - Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - element_maxvel = element->inputs->Max(VelEnum); - if(element_maxvel>maxvel) maxvel=element_maxvel; - } - - /*Figure out maximum across the cluster: */ - ISSM_MPI_Reduce(&maxvel,&node_maxvel,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&node_maxvel,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - maxvel=node_maxvel; - - /*Assign output pointers:*/ - *pmaxvel=maxvel; - -}/*}}}*/ -void FemModel::MaxVxx(IssmDouble* pmaxvx){/*{{{*/ - - int i; - IssmDouble maxvx; - IssmDouble node_maxvx; - IssmDouble element_maxvx; - - /*Go through elements, and request velocity: */ - maxvx=-INFINITY; - for(i=0;ielements->Size();i++){ - Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - element_maxvx = element->inputs->Max(VxEnum); - if(element_maxvx>maxvx) maxvx=element_maxvx; - } - - /*Figure out maximum across the cluster: */ - ISSM_MPI_Reduce(&maxvx,&node_maxvx,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&node_maxvx,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - maxvx=node_maxvx; - - /*Assign output pointers:*/ - *pmaxvx=maxvx; - -}/*}}}*/ -void FemModel::MaxVyx(IssmDouble* pmaxvy){/*{{{*/ - - int i; - IssmDouble maxvy; - IssmDouble node_maxvy; - IssmDouble element_maxvy; - - /*Go through elements, and request velocity: */ - maxvy=-INFINITY; - for(i=0;ielements->Size();i++){ - Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - element_maxvy = element->inputs->Max(VyEnum); - if(element_maxvy>maxvy) maxvy=element_maxvy; - } - - /*Figure out maximum across the cluster: */ - ISSM_MPI_Reduce(&maxvy,&node_maxvy,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&node_maxvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - maxvy=node_maxvy; - - /*Assign output pointers:*/ - *pmaxvy=maxvy; - -}/*}}}*/ -void FemModel::MaxVzx(IssmDouble* pmaxvz){/*{{{*/ - - int i; - IssmDouble maxvz; - IssmDouble node_maxvz; - IssmDouble element_maxvz; - - /*Go through elements, and request velocity: */ - maxvz=-INFINITY; - for(i=0;ielements->Size();i++){ - Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - element_maxvz = element->inputs->Max(VzEnum); - if(element_maxvz>maxvz) maxvz=element_maxvz; - } - - /*Figure out maximum across the cluster: */ - ISSM_MPI_Reduce(&maxvz,&node_maxvz,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&node_maxvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - maxvz=node_maxvz; - - /*Assign output pointers:*/ - *pmaxvz=maxvz; - -}/*}}}*/ -void FemModel::MinVelx(IssmDouble* pminvel){/*{{{*/ - - int i; - IssmDouble minvel; - IssmDouble node_minvel; - IssmDouble element_minvel; - - /*Go through elements, and request velocity: */ - minvel=INFINITY; - for(i=0;ielements->Size();i++){ - Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - element_minvel = element->inputs->Min(VelEnum); - if(element_minvelelements->Size();i++){ - Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - element_minvx = element->inputs->Min(VxEnum); - if(element_minvxelements->Size();i++){ - Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - element_minvy = element->inputs->Min(VyEnum); - if(element_minvyelements->Size();i++){ - Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - element_minvz = element->inputs->Min(VzEnum); - if(element_minvzelements->Size();i++){ - Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - local_fbmb+=element->TotalFloatingBmb(); - } - ISSM_MPI_Reduce(&local_fbmb,&total_fbmb,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&total_fbmb,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - - /*Assign output pointers: */ - *pFbmb=total_fbmb; - -}/*}}}*/ -void FemModel::TotalGroundedBmbx(IssmDouble* pGbmb){/*{{{*/ - - IssmDouble local_gbmb = 0; - IssmDouble total_gbmb; - - for(int i=0;ielements->Size();i++){ - Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - local_gbmb+=element->TotalGroundedBmb(); - } - ISSM_MPI_Reduce(&local_gbmb,&total_gbmb,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&total_gbmb,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - - /*Assign output pointers: */ - *pGbmb=total_gbmb; - -}/*}}}*/ -void FemModel::TotalSmbx(IssmDouble* pSmb){/*{{{*/ - - IssmDouble local_smb = 0; - IssmDouble total_smb; - - for(int i=0;ielements->Size();i++){ - Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - local_smb+=element->TotalSmb(); - } - ISSM_MPI_Reduce(&local_smb,&total_smb,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&total_smb,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - - /*Assign output pointers: */ - *pSmb=total_smb; - -}/*}}}*/ -void FemModel::Divergencex(IssmDouble* pdiv){/*{{{*/ - - IssmDouble local_divergence=0; - IssmDouble total_divergence; - - for(int i=0;ielements->Size();i++){ - Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - local_divergence+=element->Divergence(); - } - ISSM_MPI_Reduce(&local_divergence,&total_divergence,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&total_divergence,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - - /*Assign output pointers: */ - *pdiv=total_divergence; - -}/*}}}*/ -void FemModel::MaxDivergencex(IssmDouble* pdiv){/*{{{*/ - - IssmDouble local_divergence; - IssmDouble node_max_divergence; - IssmDouble max_divergence = -INFINITY; - - for(int i=0;ielements->Size();i++){ - Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - local_divergence=element->Divergence(); - if(fabs(local_divergence)>max_divergence) max_divergence=fabs(local_divergence); - } - ISSM_MPI_Reduce(&max_divergence,&node_max_divergence,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&node_max_divergence,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - max_divergence=node_max_divergence; - - /*Assign output pointers: */ - *pdiv=max_divergence; - -}/*}}}*/ -void FemModel::FloatingAreax(IssmDouble* pV){/*{{{*/ - - IssmDouble local_floating_area= 0; - IssmDouble total_floating_area; - - for(int i=0;ielements->Size();i++){ - Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - local_floating_area+=element->FloatingArea(); - } - ISSM_MPI_Reduce(&local_floating_area,&total_floating_area,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&total_floating_area,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - - /*Assign output pointers: */ - *pV=total_floating_area; - -}/*}}}*/ -void FemModel::GetInputLocalMinMaxOnNodesx(IssmDouble** pmin,IssmDouble** pmax,IssmDouble* ug){/*{{{*/ - - /*Get vector sizes for current configuration*/ - int configuration_type; - this->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); - int numnodes = this->nodes->NumberOfNodes(configuration_type); - - /*Initialize output vectors*/ - IssmDouble* uLmin_local = xNew(numnodes); - IssmDouble* uLmax_local = xNew(numnodes); - IssmDouble* uLmin = xNew(numnodes); - IssmDouble* uLmax = xNew(numnodes); - for(int i=0;ielements->Size();i++){ - Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - element->GetInputLocalMinMaxOnNodes(uLmin_local,uLmax_local,ug); - } - - /*Synchronize all CPUs*/ - ISSM_MPI_Allreduce((void*)uLmin_local,(void*)uLmin,numnodes,ISSM_MPI_DOUBLE,ISSM_MPI_MIN,IssmComm::GetComm()); - ISSM_MPI_Allreduce((void*)uLmax_local,(void*)uLmax,numnodes,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,IssmComm::GetComm()); - xDelete(uLmin_local); - xDelete(uLmax_local); - - /*Assign output pointers: */ - *pmin=uLmin; - *pmax=uLmax; - -}/*}}}*/ -void FemModel::GroundedAreax(IssmDouble* pV){/*{{{*/ - - IssmDouble local_grounded_area= 0; - IssmDouble total_grounded_area; - - for(int i=0;ielements->Size();i++){ - Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - local_grounded_area+=element->GroundedArea(); - } - ISSM_MPI_Reduce(&local_grounded_area,&total_grounded_area,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&total_grounded_area,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - - /*Assign output pointers: */ - *pV=total_grounded_area; - -}/*}}}*/ -void FemModel::IceMassx(IssmDouble* pM){/*{{{*/ - - IssmDouble local_ice_mass = 0; - IssmDouble total_ice_mass; - - for(int i=0;ielements->Size();i++){ - Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - local_ice_mass+=element->IceMass(); - } - ISSM_MPI_Reduce(&local_ice_mass,&total_ice_mass,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&total_ice_mass,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - - /*Assign output pointers: */ - *pM=total_ice_mass; - -}/*}}}*/ -void FemModel::IceVolumex(IssmDouble* pV){/*{{{*/ - - IssmDouble local_ice_volume = 0; - IssmDouble total_ice_volume; - - for(int i=0;ielements->Size();i++){ - Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - local_ice_volume+=element->IceVolume(); - } - ISSM_MPI_Reduce(&local_ice_volume,&total_ice_volume,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&total_ice_volume,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - - /*Assign output pointers: */ - *pV=total_ice_volume; - -}/*}}}*/ -void FemModel::IceVolumeAboveFloatationx(IssmDouble* pV){/*{{{*/ - - IssmDouble local_ice_volume_af = 0; - IssmDouble total_ice_volume_af; - - for(int i=0;ielements->Size();i++){ - Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - local_ice_volume_af+=element->IceVolumeAboveFloatation(); - } - ISSM_MPI_Reduce(&local_ice_volume_af,&total_ice_volume_af,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&total_ice_volume_af,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - - /*Assign output pointers: */ - *pV=total_ice_volume_af; - -}/*}}}*/ -void FemModel::ElementResponsex(IssmDouble* presponse,int response_enum){/*{{{*/ - - int found=0; - int sumfound=0; - int cpu_found=-1; - int index; - IssmDouble response; - Element* element=NULL; - - /*retrieve element we are interested in: */ - this->parameters->FindParam(&index,IndexEnum); - int my_rank=IssmComm::GetRank(); - - /*now, go through our elements, and retrieve the one with this id: index: */ - for(int i=0;ielements->Size();i++){ - element=xDynamicCast(this->elements->GetObjectByOffset(i)); - if (element->Id()==index){ - found=1; - cpu_found=my_rank; - break; - } - } - - /*Broadcast whether we found the element: */ - ISSM_MPI_Allreduce ( &found,&sumfound,1,ISSM_MPI_INT,ISSM_MPI_SUM,IssmComm::GetComm()); - if(!sumfound)_error_("could not find material with id" << index << " to compute ElementResponse"); - - /*Ok, we found the element, compute responseocity: */ - if(my_rank==cpu_found){ - element->ElementResponse(&response,response_enum); - } - - /*Broadcast and plug into response: */ - ISSM_MPI_Allreduce ( &cpu_found,&cpu_found,1,ISSM_MPI_INT,ISSM_MPI_MAX,IssmComm::GetComm()); - ISSM_MPI_Bcast(&response,1,ISSM_MPI_DOUBLE,cpu_found,IssmComm::GetComm()); - - /*Assign output pointers: */ - *presponse=response; - -}/*}}}*/ -void FemModel::BalancethicknessMisfitx(IssmDouble* presponse){/*{{{*/ - - /*output: */ - IssmDouble J=0.; - IssmDouble J_sum; - - IssmDouble weight,vx,vy,H,dvx[2],dvy[2],dH[2]; - IssmDouble temp,Jdet,dhdt,groundedice_melting,surface_mass_balance; - IssmDouble* xyz_list = NULL; - IssmDouble dp[3]; - - /*Compute Misfit: */ + Element* element=dynamic_cast(this->elements->GetObjectByOffset(i)); + element->StrainRateperpendicular(); + } +} +/*}}}*/ +void FemModel::StressIntensityFactorx(){/*{{{*/ + + /*Update input for basal element only*/ for(int i=0;iSize();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - - /*If on water, return 0: */ - if(!element->IsIceInElement()) continue; - - /* Get node coordinates*/ - element->GetVerticesCoordinates(&xyz_list); - Input* weights_input = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); - Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); - Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); - Input* surface_mass_balance_input = element->GetInput(SmbMassBalanceEnum); _assert_(surface_mass_balance_input); - Input* groundedice_melting_input = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(groundedice_melting_input); - Input* dhdt_input = element->GetInput(BalancethicknessThickeningRateEnum); _assert_(dhdt_input); - - /* Start looping on the number of gaussian points: */ - Gauss* gauss=element->NewGauss(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - /* Get Jacobian determinant: */ - element->JacobianDeterminant(&Jdet,xyz_list,gauss); - - /*Get all parameters at gaussian point*/ - weights_input->GetInputValue(&weight,gauss,BalancethicknessMisfitEnum); - thickness_input->GetInputValue(&H, gauss); - thickness_input->GetInputDerivativeValue(&dH[0],xyz_list,gauss); - surface_mass_balance_input->GetInputValue(&surface_mass_balance,gauss); - groundedice_melting_input->GetInputValue(&groundedice_melting,gauss); - dhdt_input->GetInputValue(&dhdt,gauss); - vx_input->GetInputValue(&vx,gauss); - vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); - vy_input->GetInputValue(&vy,gauss); - vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); - - /*Balance thickness soft constraint J = 1/2 (div(Hv)-a)^2*/ - temp = vx*dH[0]+vy*dH[1]+H*(dvx[0]+dvy[1]) - (surface_mass_balance-groundedice_melting-dhdt); - J +=weight*1/2*temp*temp*Jdet*gauss->weight; - } - - /*clean up and Return: */ - xDelete(xyz_list); - delete gauss; - } - - /*Sum all J from all cpus of the cluster:*/ - ISSM_MPI_Reduce (&J,&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - J=J_sum; - - /*Assign output pointers: */ - *presponse=J; - -}/*}}}*/ + Element* element=dynamic_cast(this->elements->GetObjectByOffset(i)); + element->StressIntensityFactor(); + } +} + /*}}}*/ void FemModel::SurfaceAbsMisfitx(IssmDouble* presponse){/*{{{*/ @@ -2059,133 +2029,163 @@ } /*}}}*/ -void FemModel::CostFunctionx(IssmDouble* pJ,IssmDouble** pJlist,int* pn){/*{{{*/ - - /*Intermediary*/ - int num_responses; - int *responses = NULL; - Results *cost_functions = NULL; - - /*Recover parameters*/ - parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum); - parameters->FindParam(&responses,NULL,InversionCostFunctionsEnum); - - /*Get the value of all cost functions*/ - this->RequestedOutputsx(&cost_functions,responses,num_responses); - - /*Get and add all contributions one by one*/ - IssmDouble J=0; - IssmDouble* Jlist = xNew(num_responses); - for(int i=0;iGetObjectByOffset(i); - Jlist[i] = reCast(result->GetValue()); - J += Jlist[i]; - } - _assert_(cost_functions->Size()==num_responses); +void FemModel::TimeAdaptx(IssmDouble* pdt){/*{{{*/ + + int i; + + /*output: */ + IssmDouble dt; + + /*intermediary: */ + Element *element = NULL; + IssmDouble min_dt = 0; + IssmDouble node_min_dt = 0; + + /*Go through elements, and figure out the minimum of the time steps for each element (using CFL criterion): */ + element=(Element*)elements->GetObjectByOffset(0); min_dt=element->TimeAdapt(); + + for (i=1;iSize();i++){ + element=xDynamicCast(elements->GetObjectByOffset(i)); + dt=element->TimeAdapt(); + if(dtelements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + local_fbmb+=element->TotalFloatingBmb(); + } + ISSM_MPI_Reduce(&local_fbmb,&total_fbmb,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&total_fbmb,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); /*Assign output pointers: */ - delete cost_functions; - xDelete(responses); - if(pJ) *pJ = J; - if(pJlist) *pJlist = Jlist; - else xDelete(Jlist); - if(pn) *pn = num_responses; -} -/*}}}*/ -void FemModel::OutputControlsx(Results **presults){/*{{{*/ - - /*parameters: */ - int num_controls,step; - IssmDouble time; - int *control_type = NULL; - - /*recover results*/ - Results* results = *presults; - if(!results) results = new Results(); - - /*Get list of Controls*/ - this->parameters->FindParam(&num_controls,InversionNumControlParametersEnum); - this->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum); - this->parameters->FindParam(&step,StepEnum); - this->parameters->FindParam(&time,TimeEnum); - - for(int i=0;i *vector_control = new Vector(this->vertices->NumberOfVertices()); - Vector *vector_gradient = new Vector(this->vertices->NumberOfVertices()); - - /*Fill in vector*/ - for(int j=0;jSize();j++){ - Element* element=(Element*)elements->GetObjectByOffset(j); - element->ControlToVectors(vector_control,vector_gradient,control_enum); - } - vector_control->Assemble(); - vector_gradient->Assemble(); - - results->AddResult(new GenericExternalResult*>(results->Size()+1,control_enum,vector_control ,step,time)); - results->AddResult(new GenericExternalResult*>(results->Size()+1,gradient_enum,vector_gradient,step,time)); - } - - /*Clean up and return*/ - xDelete(control_type); -} -/*}}}*/ -void FemModel::StressIntensityFactorx(){/*{{{*/ - - /*Update input for basal element only*/ + *pFbmb=total_fbmb; + +}/*}}}*/ +void FemModel::TotalGroundedBmbx(IssmDouble* pGbmb){/*{{{*/ + + IssmDouble local_gbmb = 0; + IssmDouble total_gbmb; + + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + local_gbmb+=element->TotalGroundedBmb(); + } + ISSM_MPI_Reduce(&local_gbmb,&total_gbmb,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&total_gbmb,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + + /*Assign output pointers: */ + *pGbmb=total_gbmb; + +}/*}}}*/ +void FemModel::TotalSmbx(IssmDouble* pSmb){/*{{{*/ + + IssmDouble local_smb = 0; + IssmDouble total_smb; + + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + local_smb+=element->TotalSmb(); + } + ISSM_MPI_Reduce(&local_smb,&total_smb,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&total_smb,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + + /*Assign output pointers: */ + *pSmb=total_smb; + +}/*}}}*/ +void FemModel::UpdateConstraintsExtrudeFromBasex(void){ /*{{{*/ + for(int i=0;iSize();i++){ - Element* element=dynamic_cast(this->elements->GetObjectByOffset(i)); - element->StressIntensityFactor(); - } -} - /*}}}*/ -void FemModel::CalvingRateLevermannx(){/*{{{*/ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + element->UpdateConstraintsExtrudeFromBase(); + } + +} +/*}}}*/ +void FemModel::UpdateConstraintsExtrudeFromTopx(void){ /*{{{*/ for(int i=0;iSize();i++){ - Element* element=dynamic_cast(this->elements->GetObjectByOffset(i)); - element->CalvingRateLevermann(); - } -} -/*}}}*/ -void FemModel::CalvingRateDevx(){/*{{{*/ - - for(int i=0;iSize();i++){ - Element* element=dynamic_cast(this->elements->GetObjectByOffset(i)); - element->CalvingRateDev(); - } -} -/*}}}*/ -void FemModel::StrainRateparallelx(){/*{{{*/ - - for(int i=0;iSize();i++){ - Element* element=dynamic_cast(this->elements->GetObjectByOffset(i)); - element->StrainRateparallel(); - } -} -/*}}}*/ -void FemModel::StrainRateperpendicularx(){/*{{{*/ - - for(int i=0;iSize();i++){ - Element* element=dynamic_cast(this->elements->GetObjectByOffset(i)); - element->StrainRateperpendicular(); - } -} -/*}}}*/ -void FemModel::DeviatoricStressx(){/*{{{*/ - - for(int i=0;iSize();i++){ - Element* element=dynamic_cast(this->elements->GetObjectByOffset(i)); - element->ComputeDeviatoricStressTensor(); - } + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + element->UpdateConstraintsExtrudeFromTop(); + } + +} +/*}}}*/ +void FemModel::UpdateConstraintsx(void){ /*{{{*/ + + IssmDouble time; + int analysis_type; + + /*retrieve parameters: */ + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + parameters->FindParam(&time,TimeEnum); + + /*start module: */ + if(VerboseModule()) _printf0_(" Updating constraints and active domain of analysis " << EnumToStringx(analysis_type) << " for time: " << time << "\n"); + + Analysis* analysis= EnumToAnalysis(analysis_type); + analysis->UpdateConstraints(this); + delete analysis; + + /*Second, constraints might be time dependent: */ + SpcNodesx(nodes,constraints,parameters,analysis_type); + + /*Now, update degrees of freedoms: */ + NodesDofx(nodes,parameters,analysis_type); + +} +/*}}}*/ +int FemModel::UpdateVertexPositionsx(void){ /*{{{*/ + + IssmDouble *surface = NULL; + IssmDouble *bed = NULL; + + if(VerboseSolution()) _printf0_(" updating vertices positions\n"); + + /*get vertex vectors for bed and thickness: */ + GetVectorFromInputsx(&surface ,this, SurfaceEnum,VertexPIdEnum); + GetVectorFromInputsx(&bed ,this, BaseEnum, VertexPIdEnum); + + /*Allocate vector*/ + Vector *vx=new Vector(vertices->NumberOfVertices()); + Vector *vy=new Vector(vertices->NumberOfVertices()); + Vector *vz=new Vector(vertices->NumberOfVertices()); + + /*Update verices new geometry: */ + for(int i=0;iSize();i++){ + Vertex* vertex=(Vertex*)vertices->GetObjectByOffset(i); + vertex->UpdatePosition(vx,vy,vz,parameters,surface,bed); + } + + /*Assemble mesh velocity*/ + vx->Assemble(); + vy->Assemble(); + vz->Assemble(); + + /*Update element inputs*/ + InputUpdateFromVectorx(this,vx,VxMeshEnum,VertexPIdEnum); + InputUpdateFromVectorx(this,vy,VyMeshEnum,VertexPIdEnum); + InputUpdateFromVectorx(this,vz,VzMeshEnum,VertexPIdEnum); + + /*Free ressources:*/ + delete vx; + delete vy; + delete vz; + xDelete(bed); + xDelete(surface); + return 1; } /*}}}*/ @@ -2608,74 +2608,4 @@ } /*}}}*/ -void FemModel::InitFromBuffers(char* buffer, int buffersize, char* toolkits, int in_solution_type, bool trace, IssmPDouble* X){/*{{{*/ - - /*intermediary*/ - int i; - int analysis_type; - FILE *IOMODEL = NULL; - FILE *toolkitsoptionsfid = NULL; - FILE *output_fid = NULL; - int my_rank; - size_t outputsize; - char *outputbuffer; - const char* rootpath=""; //needed for Dakota runs only, which we won't do here. - - /*recover my_rank:*/ - my_rank=IssmComm::GetRank(); - - /*Initialize internal data: */ - this->solution_type = in_solution_type; - this->analysis_counter = nummodels-1; //point to last analysis_type carried out. - this->results = new Results(); //not initialized by CreateDataSets - - /*Open input file on cpu 0 and create IoModel */ - if(my_rank==0) IOMODEL = fmemopen((void*)buffer, buffersize, "rb"); - IoModel* iomodel = new IoModel(IOMODEL,in_solution_type,trace,X); - - /*Figure out what analyses are activated for this solution*/ - SolutionAnalysesList(&this->analysis_type_list,&this->nummodels,iomodel,this->solution_type); - - /*Open toolkits file: */ - toolkitsoptionsfid=fmemopen((void*)toolkits, strlen(toolkits)+1, "r"); - - /*create datasets for all analyses*/ - ModelProcessorx(&this->elements,&this->nodes,&this->vertices,&this->materials,&this->constraints,&this->loads,&this->parameters,iomodel,toolkitsoptionsfid,(char*)rootpath,this->solution_type,this->nummodels,this->analysis_type_list); - - /*do the post-processing of the datasets to get an FemModel that can actually run analyses: */ - for(i=0;iSetCurrentConfiguration(analysis_type); - - if(i==0){ - if(VerboseMProcessor()) _printf0_(" creating vertex PIDs\n"); - VerticesDofx(vertices,parameters); //only call once, we only have one set of vertices - } - - if(VerboseMProcessor()) _printf0_(" resolving node constraints\n"); - SpcNodesx(nodes,constraints,parameters,analysis_type); - - if(VerboseMProcessor()) _printf0_(" creating nodal degrees of freedom\n"); - NodesDofx(nodes,parameters,analysis_type); - - if(VerboseMProcessor()) _printf0_(" configuring element and loads\n"); - ConfigureObjectsx(elements, loads, nodes, vertices, materials,parameters); - } - - /*Close input file and toolkits file descriptors: */ - if(my_rank==0) fclose(IOMODEL); - fclose(toolkitsoptionsfid); - - /*Open output file once for all and add output file descriptor to parameters*/ - output_fid=open_memstream(&outputbuffer,&outputsize); - if(output_fid==NULL)_error_("could not initialize output stream"); - this->parameters->SetParam(output_fid,OutputFilePointerEnum); - this->parameters->AddObject(new GenericParam(&outputbuffer,OutputBufferPointerEnum)); - this->parameters->AddObject(new GenericParam(&outputsize,OutputBufferSizePointerEnum)); - - /*Clean up*/ - delete iomodel; -}/*}}}*/ void FemModel::CleanUpJs(char** poutput, size_t* psize){/*{{{*/ @@ -2717,3 +2647,73 @@ } /*}}}*/ +void FemModel::InitFromBuffers(char* buffer, int buffersize, char* toolkits, int in_solution_type, bool trace, IssmPDouble* X){/*{{{*/ + + /*intermediary*/ + int i; + int analysis_type; + FILE *IOMODEL = NULL; + FILE *toolkitsoptionsfid = NULL; + FILE *output_fid = NULL; + int my_rank; + size_t outputsize; + char *outputbuffer; + const char* rootpath=""; //needed for Dakota runs only, which we won't do here. + + /*recover my_rank:*/ + my_rank=IssmComm::GetRank(); + + /*Initialize internal data: */ + this->solution_type = in_solution_type; + this->analysis_counter = nummodels-1; //point to last analysis_type carried out. + this->results = new Results(); //not initialized by CreateDataSets + + /*Open input file on cpu 0 and create IoModel */ + if(my_rank==0) IOMODEL = fmemopen((void*)buffer, buffersize, "rb"); + IoModel* iomodel = new IoModel(IOMODEL,in_solution_type,trace,X); + + /*Figure out what analyses are activated for this solution*/ + SolutionAnalysesList(&this->analysis_type_list,&this->nummodels,iomodel,this->solution_type); + + /*Open toolkits file: */ + toolkitsoptionsfid=fmemopen((void*)toolkits, strlen(toolkits)+1, "r"); + + /*create datasets for all analyses*/ + ModelProcessorx(&this->elements,&this->nodes,&this->vertices,&this->materials,&this->constraints,&this->loads,&this->parameters,iomodel,toolkitsoptionsfid,(char*)rootpath,this->solution_type,this->nummodels,this->analysis_type_list); + + /*do the post-processing of the datasets to get an FemModel that can actually run analyses: */ + for(i=0;iSetCurrentConfiguration(analysis_type); + + if(i==0){ + if(VerboseMProcessor()) _printf0_(" creating vertex PIDs\n"); + VerticesDofx(vertices,parameters); //only call once, we only have one set of vertices + } + + if(VerboseMProcessor()) _printf0_(" resolving node constraints\n"); + SpcNodesx(nodes,constraints,parameters,analysis_type); + + if(VerboseMProcessor()) _printf0_(" creating nodal degrees of freedom\n"); + NodesDofx(nodes,parameters,analysis_type); + + if(VerboseMProcessor()) _printf0_(" configuring element and loads\n"); + ConfigureObjectsx(elements, loads, nodes, vertices, materials,parameters); + } + + /*Close input file and toolkits file descriptors: */ + if(my_rank==0) fclose(IOMODEL); + fclose(toolkitsoptionsfid); + + /*Open output file once for all and add output file descriptor to parameters*/ + output_fid=open_memstream(&outputbuffer,&outputsize); + if(output_fid==NULL)_error_("could not initialize output stream"); + this->parameters->SetParam(output_fid,OutputFilePointerEnum); + this->parameters->AddObject(new GenericParam(&outputbuffer,OutputBufferPointerEnum)); + this->parameters->AddObject(new GenericParam(&outputsize,OutputBufferSizePointerEnum)); + + /*Clean up*/ + delete iomodel; +}/*}}}*/ #endif Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 20810) @@ -28,19 +28,19 @@ public: + int analysis_counter; //counter into analysis_type_list + int *analysis_type_list; //list of analyses this femmodel is going to carry out int nummodels; int solution_type; - int *analysis_type_list; //list of analyses this femmodel is going to carry out - int analysis_counter; //counter into analysis_type_list Profiler* profiler; //keep time, cpu and mem statistics while we are running. + Constraints *constraints; //one set of constraints. each constraint knows which analysis_type it handles Elements *elements; //elements (one set for all analyses) - Nodes *nodes; //one set of nodes - Vertices *vertices; //one set of vertices - Constraints *constraints; //one set of constraints. each constraint knows which analysis_type it handles Loads *loads; //one set of constraints. each constraint knows which analysis_type it handles Materials *materials; //one set of materials, for each element + Nodes *nodes; //one set of nodes Parameters *parameters; //one set of parameters, independent of the analysis_type Results *results; //results that cannot be fit into the elements + Vertices *vertices; //one set of vertices /*constructors, destructors: */ @@ -50,23 +50,35 @@ /*Methods:*/ + void CheckPoint(void); + void CleanUp(void); + FemModel* copy(); void Echo(); - FemModel* copy(); void InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, char* restartfilename, const int solution_type,bool trace,IssmPDouble* X=NULL); - void SolutionAnalysesList(int** panalyses,int* pnumanalyses,IoModel* iomodel,int solutiontype); - void CleanUp(void); - void Solve(void); + void Marshall(char** pmarshalled_data, int* pmarshalled_data_size, int marshall_direction); + void Restart(void); void SetCurrentConfiguration(int configuration_type); void SetCurrentConfiguration(int configuration_type,int analysis_type); - void CheckPoint(void); - void Marshall(char** pmarshalled_data, int* pmarshalled_data_size, int marshall_direction); - void Restart(void); + void SolutionAnalysesList(int** panalyses,int* pnumanalyses,IoModel* iomodel,int solutiontype); + void Solve(void); /*Modules*/ + void BalancethicknessMisfitx(IssmDouble* pV); + void CalvingRateDevx(); + void CalvingRateLevermannx(); + void DeviatoricStressx(); + void Divergencex(IssmDouble* pdiv); void ElementOperationx(void (Element::*function)(void)); + void ElementResponsex(IssmDouble* presponse,int response_enum); + void FloatingAreax(IssmDouble* pV); void GetInputLocalMinMaxOnNodesx(IssmDouble** pmin,IssmDouble** pmax,IssmDouble* ug); + void GroundedAreax(IssmDouble* pV); + void IceMassx(IssmDouble* pV); + void IceVolumex(IssmDouble* pV); + void IceVolumeAboveFloatationx(IssmDouble* pV); void MassFluxx(IssmDouble* presponse); void MaxAbsVxx(IssmDouble* presponse); void MaxAbsVyx(IssmDouble* presponse); void MaxAbsVzx(IssmDouble* presponse); + void MaxDivergencex(IssmDouble* pdiv); void MaxVelx(IssmDouble* presponse); void MaxVxx(IssmDouble* presponse); @@ -77,36 +89,24 @@ void MinVyx(IssmDouble* presponse); void MinVzx(IssmDouble* presponse); + void ResetLevelset(); + void StrainRateparallelx(); + void StrainRateperpendicularx(); + void StressIntensityFactorx(); void TotalFloatingBmbx(IssmDouble* pFbmb); void TotalGroundedBmbx(IssmDouble* pGbmb); void TotalSmbx(IssmDouble* pSmb); - void Divergencex(IssmDouble* pdiv); - void MaxDivergencex(IssmDouble* pdiv); - void FloatingAreax(IssmDouble* pV); - void GroundedAreax(IssmDouble* pV); - void IceMassx(IssmDouble* pV); - void IceVolumex(IssmDouble* pV); - void IceVolumeAboveFloatationx(IssmDouble* pV); - void ElementResponsex(IssmDouble* presponse,int response_enum); - void BalancethicknessMisfitx(IssmDouble* pV); - void StressIntensityFactorx(); - void StrainRateparallelx(); - void StrainRateperpendicularx(); - void DeviatoricStressx(); - void CalvingRateLevermannx(); - void CalvingRateDevx(); - void ResetLevelset(); #ifdef _HAVE_DAKOTA_ void DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses); #endif + void CostFunctionx(IssmDouble* pJ,IssmDouble** pJlist,int* pn); + void OutputControlsx(Results **presults); + void RequestedDependentsx(void); void RequestedOutputsx(Results **presults,char** requested_outputs, int numoutputs,bool save_results=true); void RequestedOutputsx(Results **presults,int* requested_outputs, int numoutputs,bool save_results=true); - void RequestedDependentsx(void); void Responsex(IssmDouble* presponse,int response_descriptor_enum); void Responsex(IssmDouble* presponse,const char* response_descriptor); - void OutputControlsx(Results **presults); - void CostFunctionx(IssmDouble* pJ,IssmDouble** pJlist,int* pn); + void SurfaceAbsMisfitx( IssmDouble* pJ); void ThicknessAbsGradientx( IssmDouble* pJ); void ThicknessPositivex(IssmDouble* pJ); - void SurfaceAbsMisfitx( IssmDouble* pJ); #ifdef _HAVE_GIA_ void Deflection(Vector* wg,Vector* dwgdt, IssmDouble* x, IssmDouble* y); @@ -117,16 +117,16 @@ IssmDouble SealevelriseOceanAverage(Vector* Sg); #endif + void HydrologyEPLupdateDomainx(IssmDouble* pEplcount); void TimeAdaptx(IssmDouble* pdt); + void UpdateConstraintsExtrudeFromBasex(); + void UpdateConstraintsExtrudeFromTopx(); + void UpdateConstraintsL2ProjectionEPLx(IssmDouble* pL2count); void UpdateConstraintsx(void); int UpdateVertexPositionsx(void); - void UpdateConstraintsExtrudeFromBasex(); - void UpdateConstraintsExtrudeFromTopx(); - void HydrologyEPLupdateDomainx(IssmDouble* pEplcount); - void UpdateConstraintsL2ProjectionEPLx(IssmDouble* pL2count); #ifdef _HAVE_JAVASCRIPT_ FemModel(IssmDouble* buffer, int buffersize, char* toolkits, char* solution, char* modelname,ISSM_MPI_Comm incomm, bool trace=false); + void CleanUpJs(char** poutput, size_t* psize); void InitFromBuffers(char* buffer, int buffersize, char* toolkits, int solution_type,bool trace,IssmPDouble* X=NULL); - void CleanUpJs(char** poutput, size_t* psize); #endif }; Index: /issm/trunk-jpl/src/c/classes/GiaDeflectionCoreArgs.h =================================================================== --- /issm/trunk-jpl/src/c/classes/GiaDeflectionCoreArgs.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/GiaDeflectionCoreArgs.h (revision 20810) @@ -9,18 +9,18 @@ /*inputs: */ + IssmDouble currenttime; + IssmDouble* hes; //loading history (in ice thickness) + int numtimes; //loading history length IssmDouble ri; //radial distance from center of disk to vertex i IssmDouble re; //radius of disk - IssmDouble* hes; //loading history (in ice thickness) IssmDouble* times; //loading history times - int numtimes; //loading history length - IssmDouble currenttime; /*gia material parameters: */ + IssmDouble lithosphere_density; IssmDouble lithosphere_shear_modulus; - IssmDouble lithosphere_density; + IssmDouble lithosphere_thickness; + IssmDouble mantle_density; IssmDouble mantle_shear_modulus; IssmDouble mantle_viscosity; - IssmDouble mantle_density; - IssmDouble lithosphere_thickness; /*gia solution parameters: */ Index: /issm/trunk-jpl/src/c/classes/Hook.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Hook.cpp (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Hook.cpp (revision 20810) @@ -59,20 +59,18 @@ /*Some of the Object functionality: */ -void Hook::Echo(void){/*{{{*/ - _assert_(this); - int i; - if (num){ - _printf_(" Hook: \n"); - _printf_(" num=" << this->num << "\n"); - _printf_(" ids: "); - for(i=0;inum;i++) _printf_(this->ids[i] << " "); - _printf_("\n"); - _printf_(" offsets: "); - for (i=0;inum;i++) _printf_(this->offsets[i] << " "); - _printf_("\n"); - } - else{ - _printf_(" Hook: num=0 \n"); - } +Object* Hook::copy(void){/*{{{*/ + + /*output: */ + Hook* output=NULL; + + /*initalize output: */ + output=new Hook(this->ids,this->num); + + for(int i=0;inum;i++){ + output->objects[i] = this->objects[i]; + output->offsets[i] = this->offsets[i]; + } + + return (Object*)output; } /*}}}*/ @@ -104,18 +102,20 @@ } /*}}}*/ -Object* Hook::copy(void){/*{{{*/ - - /*output: */ - Hook* output=NULL; - - /*initalize output: */ - output=new Hook(this->ids,this->num); - - for(int i=0;inum;i++){ - output->objects[i] = this->objects[i]; - output->offsets[i] = this->offsets[i]; - } - - return (Object*)output; +void Hook::Echo(void){/*{{{*/ + _assert_(this); + int i; + if (num){ + _printf_(" Hook: \n"); + _printf_(" num=" << this->num << "\n"); + _printf_(" ids: "); + for(i=0;inum;i++) _printf_(this->ids[i] << " "); + _printf_("\n"); + _printf_(" offsets: "); + for (i=0;inum;i++) _printf_(this->offsets[i] << " "); + _printf_("\n"); + } + else{ + _printf_(" Hook: num=0 \n"); + } } /*}}}*/ @@ -143,15 +143,4 @@ /*Hook management: */ -void Hook::reset(){/*{{{*/ - - /*intermediary: */ - Object* object=NULL; - int i; - - for(i=0;inum;i++){ - this->objects[i]=NULL; //reset this node. - } -} -/*}}}*/ void Hook::configure(DataSet* dataset){/*{{{*/ @@ -196,4 +185,8 @@ } /*}}}*/ +Object** Hook::deliverp(void){/*{{{*/ + return objects; +} +/*}}}*/ Object* Hook::delivers(void){/*{{{*/ @@ -208,6 +201,6 @@ /*}}}*/ -Object** Hook::deliverp(void){/*{{{*/ - return objects; +int Hook::GetNum(void){/*{{{*/ + return this->num; } /*}}}*/ @@ -216,6 +209,13 @@ } /*}}}*/ -int Hook::GetNum(void){/*{{{*/ - return this->num; +void Hook::reset(){/*{{{*/ + + /*intermediary: */ + Object* object=NULL; + int i; + + for(i=0;inum;i++){ + this->objects[i]=NULL; //reset this node. + } } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Hook.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Hook.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Hook.h (revision 20810) @@ -18,7 +18,7 @@ private: + int *ids; //list of object ids, to go look for them in datasets. int num; //number of objects being hooked onto Object **objects; //list of object pointers - int *ids; //list of object ids, to go look for them in datasets. int *offsets; //list of object offsets into datasets, to speed up lookup. @@ -31,17 +31,17 @@ /*}}}*/ /*Object like functionality:{{{*/ + Object* copy(void); + void DeepEcho(void); void Echo(void); - void DeepEcho(void); - Object* copy(void); void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction); /*}}}*/ /*Hook management: {{{*/ + void configure(DataSet* dataset); + Object** deliverp(void); //deliver all objects Object* delivers(void); //single object deliver - Object** deliverp(void); //deliver all objects - void configure(DataSet* dataset); + int GetNum(void); + int* Ids(void); void reset(void); Hook* Spawn(int* indices, int numindices); - int* Ids(void); - int GetNum(void); /*}}}*/ }; Index: /issm/trunk-jpl/src/c/classes/Inputs/BoolInput.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/BoolInput.cpp (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Inputs/BoolInput.cpp (revision 20810) @@ -29,6 +29,8 @@ /*Object virtual functions definitions:*/ -void BoolInput::Echo(void){/*{{{*/ - this->DeepEcho(); +Object* BoolInput::copy() {/*{{{*/ + + return new BoolInput(this->enum_type,this->value); + } /*}}}*/ @@ -38,17 +40,9 @@ } /*}}}*/ -int BoolInput::Id(void){ return -1; }/*{{{*/ -/*}}}*/ -int BoolInput::ObjectEnum(void){/*{{{*/ - - return BoolInputEnum; - +void BoolInput::Echo(void){/*{{{*/ + this->DeepEcho(); } /*}}}*/ -Object* BoolInput::copy() {/*{{{*/ - - return new BoolInput(this->enum_type,this->value); - -} +int BoolInput::Id(void){ return -1; }/*{{{*/ /*}}}*/ void BoolInput::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/ @@ -61,4 +55,10 @@ } /*}}}*/ +int BoolInput::ObjectEnum(void){/*{{{*/ + + return BoolInputEnum; + +} +/*}}}*/ /*BoolInput management*/ @@ -66,4 +66,18 @@ return this->enum_type; + +} +/*}}}*/ +Input* BoolInput::SpawnSegInput(int index1,int index2){/*{{{*/ + + /*output*/ + BoolInput* outinput=new BoolInput(); + + /*only copy current value*/ + outinput->enum_type=this->enum_type; + outinput->value=this->value; + + /*Assign output*/ + return outinput; } @@ -83,43 +97,6 @@ } /*}}}*/ -Input* BoolInput::SpawnSegInput(int index1,int index2){/*{{{*/ - - /*output*/ - BoolInput* outinput=new BoolInput(); - - /*only copy current value*/ - outinput->enum_type=this->enum_type; - outinput->value=this->value; - - /*Assign output*/ - return outinput; - -} -/*}}}*/ /*Object functions*/ -void BoolInput::GetInputValue(bool* pvalue){/*{{{*/ - *pvalue=value; -} -/*}}}*/ -void BoolInput::GetInputValue(int* pvalue){_error_("not supported yet!");}/*{{{*/ -/*}}}*/ -void BoolInput::GetInputValue(IssmDouble* pvalue){_error_("not supported yet!");}/*{{{*/ -/*}}}*/ -void BoolInput::GetInputValue(IssmDouble* pvalue,Gauss* gauss){_error_("not supported yet!");}/*{{{*/ -/*}}}*/ -void BoolInput::ChangeEnum(int newenumtype){/*{{{*/ - this->enum_type=newenumtype; -} -/*}}}*/ -void BoolInput::SquareMin(IssmDouble* psquaremin,Parameters* parameters){/*{{{*/ - /*square of a bool is the bool itself: */ - *psquaremin=value; -} -/*}}}*/ -void BoolInput::Scale(IssmDouble scale_factor){/*{{{*/ - /*a bool cannot be scaled: */ -} -/*}}}*/ void BoolInput::AXPY(Input* xinput,IssmDouble scalar){/*{{{*/ @@ -142,4 +119,12 @@ } /*}}}*/ +void BoolInput::ChangeEnum(int newenumtype){/*{{{*/ + this->enum_type=newenumtype; +} +/*}}}*/ +void BoolInput::Configure(Parameters* parameters){/*{{{*/ + /*do nothing: */ +} +/*}}}*/ void BoolInput::Extrude(int start){/*{{{*/ @@ -147,4 +132,14 @@ return; } +/*}}}*/ +void BoolInput::GetInputValue(bool* pvalue){/*{{{*/ + *pvalue=value; +} +/*}}}*/ +void BoolInput::GetInputValue(int* pvalue){_error_("not supported yet!");}/*{{{*/ +/*}}}*/ +void BoolInput::GetInputValue(IssmDouble* pvalue){_error_("not supported yet!");}/*{{{*/ +/*}}}*/ +void BoolInput::GetInputValue(IssmDouble* pvalue,Gauss* gauss){_error_("not supported yet!");}/*{{{*/ /*}}}*/ void BoolInput::GetVectorFromInputs(Vector* vector,int* doflist){/*{{{*/ @@ -154,6 +149,11 @@ } /*}}}*/ -void BoolInput::Configure(Parameters* parameters){/*{{{*/ - /*do nothing: */ +void BoolInput::Scale(IssmDouble scale_factor){/*{{{*/ + /*a bool cannot be scaled: */ } /*}}}*/ +void BoolInput::SquareMin(IssmDouble* psquaremin,Parameters* parameters){/*{{{*/ + /*square of a bool is the bool itself: */ + *psquaremin=value; +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Inputs/BoolInput.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/BoolInput.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Inputs/BoolInput.h (revision 20810) @@ -22,26 +22,36 @@ /*}}}*/ /*Object virtual functions definitions:{{{ */ + Object *copy(); + void DeepEcho(); void Echo(); - void DeepEcho(); int Id(); + void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction); int ObjectEnum(); - Object *copy(); - void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction); /*}}}*/ /*BoolInput management: {{{*/ + void AddTimeValues(IssmDouble* values,int step,IssmDouble time){_error_("not supported yet");}; + void Configure(Parameters* parameters); + int GetResultArraySize(void){return 1;}; + int GetResultInterpolation(void){return P0Enum;}; + int GetResultNumberOfNodes(void){return 1;}; int InstanceEnum(); - Input* SpawnTriaInput(int index1,int index2,int index3); - Input* SpawnSegInput(int index1,int index2); Input* PointwiseDivide(Input* inputB){_error_("not implemented yet");}; Input* PointwiseMin(Input* inputB){_error_("not implemented yet");}; Input* PointwiseMax(Input* inputB){_error_("not implemented yet");}; - int GetResultInterpolation(void){return P0Enum;}; - int GetResultNumberOfNodes(void){return 1;}; - int GetResultArraySize(void){return 1;}; + Input* SpawnSegInput(int index1,int index2); + Input* SpawnTriaInput(int index1,int index2,int index3); void ResultToPatch(IssmDouble* values,int nodesperelement,int sid){_error_("not supported yet");}; - void Configure(Parameters* parameters); - void AddTimeValues(IssmDouble* values,int step,IssmDouble time){_error_("not supported yet");}; /*}}}*/ /*numerics: {{{*/ + void AXPY(Input* xinput,IssmDouble scalar); + void ChangeEnum(int newenumtype); + void Constrain(IssmDouble cm_min, IssmDouble cm_max){_error_("Constrain not implemented for booleans");}; + void ConstrainMin(IssmDouble minimum){_error_("not implemented yet");}; + void Extrude(int start); + void GetInputAllTimeAverages(IssmDouble** pvalues,IssmDouble** ptimes, int* pnumtimes){_error_("not implemented yet");}; + void GetInputAverage(IssmDouble* pvalue){_error_("not implemented yet");}; + void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, Gauss* gauss){_error_("not implemented yet");}; + void GetInputDerivativeAverageValue(IssmDouble* derivativevalues, IssmDouble* xyz_list){_error_("not implemented yet");}; + void GetInputUpToCurrentTimeAverages(IssmDouble** pvalues, IssmDouble** ptimes, int* pnumtimes, IssmDouble currenttime){_error_("not implemented yet");}; void GetInputValue(bool* pvalue); void GetInputValue(int* pvalue); @@ -50,12 +60,5 @@ void GetInputValue(IssmDouble* pvalue,Gauss* gauss,IssmDouble time){_error_("not implemented yet");}; void GetInputValue(IssmDouble* pvalue,Gauss* gauss ,int index){_error_("not implemented yet");}; - void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, Gauss* gauss){_error_("not implemented yet");}; - void GetInputAverage(IssmDouble* pvalue){_error_("not implemented yet");}; - void GetInputDerivativeAverageValue(IssmDouble* derivativevalues, IssmDouble* xyz_list){_error_("not implemented yet");}; - void GetInputAllTimeAverages(IssmDouble** pvalues,IssmDouble** ptimes, int* pnumtimes){_error_("not implemented yet");}; - void GetInputUpToCurrentTimeAverages(IssmDouble** pvalues, IssmDouble** ptimes, int* pnumtimes, IssmDouble currenttime){_error_("not implemented yet");}; - void ChangeEnum(int newenumtype); - void SquareMin(IssmDouble* psquaremin, Parameters* parameters); - void ConstrainMin(IssmDouble minimum){_error_("not implemented yet");}; + void GetVectorFromInputs(Vector* vector,int* doflist); IssmDouble InfinityNorm(void){_error_("InfinityNorm not implemented for booleans");}; IssmDouble Max(void){_error_("Max not implemented for booleans");}; @@ -65,9 +68,6 @@ void Set(IssmDouble setvalue){_error_("Set not implemented yet");}; void Scale(IssmDouble scale_factor); - void AXPY(Input* xinput,IssmDouble scalar); - void Constrain(IssmDouble cm_min, IssmDouble cm_max){_error_("Constrain not implemented for booleans");}; - void Extrude(int start); + void SquareMin(IssmDouble* psquaremin, Parameters* parameters); void VerticallyIntegrate(Input* thickness_input){_error_("not supported yet");}; - void GetVectorFromInputs(Vector* vector,int* doflist); /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/IoModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/IoModel.cpp (revision 20809) +++ /issm/trunk-jpl/src/c/classes/IoModel.cpp (revision 20810) @@ -377,248 +377,4 @@ } /*}}}*/ -void IoModel::FetchIndependentConstant(int* pXcount,IssmPDouble* X,const char* constant_name){/*{{{*/ - - /*recover my_rank:*/ - int my_rank=IssmComm::GetRank(); - - /*recover Xcount if X is not NULL:*/ - int Xcount = 0; - if(X) Xcount=*pXcount; - - #ifdef _HAVE_ADOLC_ //cannot come here unless you are running AD mode, from DeclaredIndependents: - - /*output: */ - IssmPDouble pscalar; - IssmDouble scalar; //same as pscalar, except it's an ADOLC independent variable - int code; - - /*Set file pointer to beginning of the data: */ - fid=this->SetFilePointerToData(&code,NULL,constant_name); - if(code!=3) _error_("expecting a IssmDouble for \"" << constant_name<<"\""); - - /*We have to read a scalar from disk. First read the dimensions of the scalar, then the scalar: */ - if(my_rank==0){ - if(fread(&pscalar,sizeof(IssmPDouble),1,fid)!=1)_error_("could not read scalar "); - - /*Now, before we even broadcast this to other nodes, declare the scalar as an independent variable!. If we - *have been supplied an X vector, use it instead of what we just read: */ - if(X){ - scalar<<=X[Xcount]; - } - else{ - scalar<<=pscalar; - } - } - - ISSM_MPI_Bcast(&scalar,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - this->AddConstantIndependent(new IoConstant(scalar,constant_name)); - - /*increment offset into X vector, now that we have read 1 value:*/ - Xcount++; *pXcount=Xcount; - #endif -} -/*}}}*/ -void IoModel::FetchIndependentData(int* pXcount,IssmPDouble* X,const char* data_name){/*{{{*/ - - /*recover my_rank:*/ - int my_rank=IssmComm::GetRank(); - - /*recover Xcount if X is not NULL:*/ - int Xcount = 0; - if(X) Xcount=*pXcount; - - #ifdef _HAVE_ADOLC_ //cannot come here unless you are running AD mode, from DeclaredIndependents: - - /*Intermediaries*/ - int M,N; - IssmPDouble* buffer=NULL; //a buffer to read the data from disk - IssmDouble* matrix=NULL; //our independent variable - int code,layout; - - /*Set file pointer to beginning of the data: */ - fid=this->SetFilePointerToData(&code,&layout,data_name); - if((code!=5) && (code!=6) && (code!=7))_error_("expecting a IssmDouble, integer or boolean matrix for \"" << data_name<<"\""); - - /*We have to read a matrix from disk. First read the dimensions of the matrix, then the whole matrix: */ - /*numberofelements: */ - if(my_rank==0){ - if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows for matrix "); - } - ISSM_MPI_Bcast(&M,1,ISSM_MPI_INT,0,IssmComm::GetComm()); - - if(my_rank==0){ - if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns for matrix "); - } - ISSM_MPI_Bcast(&N,1,ISSM_MPI_INT,0,IssmComm::GetComm()); - - /*Now allocate matrix: */ - if(M*N){ - buffer=xNew(M*N); - matrix=xNew(M*N); - - /*Read matrix on node 0, then broadcast: */ - if(my_rank==0){ - if(fread(buffer,M*N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read matrix "); - - /*Now, before we even broadcast this to other nodes, declare the whole matrix as a independent variable! - If we have been supplied an X vector, use it instead of what we just read: */ - if(X){ - for(int i=0;i(buffer); - } - else _error_("cannot declare the independent variable \"" << data_name << "\" if it's empty!"); - - /*Add to data as independent*/ - this->AddDataIndependent(new IoData(matrix,code,layout,M,N,data_name)); - - /*increment offset into X vector, now that we have read M*N values:*/ - Xcount+=M*N; *pXcount=Xcount; - #endif -} -/*}}}*/ -void IoModel::FindConstant(bool* pvalue,const char* constant_name){/*{{{*/ - - /*Intermediary*/ - vector::iterator iter; - - for(iter=constants.begin();itername,constant_name)==0){ - ioconstant->constant->GetParameterValue(pvalue); - return; - } - } - - for(vector::iterator iter=constants.begin();iterconstant->Echo(); - _error_("Could not find constant \""<::iterator iter; - - for(iter=constants.begin();itername,constant_name)==0){ - ioconstant->constant->GetParameterValue(pvalue); - return; - } - } - - _error_("Could not find constant \""<::iterator iter; - - for(iter=constants.begin();itername,constant_name)==0){ - ioconstant->constant->GetParameterValue(pvalue); - return; - } - } - - _error_("Could not find constant \""<::iterator iter; - - for(iter=constants.begin();itername,constant_name)==0){ - ioconstant->constant->GetParameterValue(pvalue); - return; - } - } - - _error_("Could not find constant \""<::iterator iter; - - for(iter=constants.begin();itername,constant_name)==0){ - ioconstant->constant->GetParameterValue(pvalue,psize); - return; - } - } - - _error_("Could not find constant \""<::iterator iter=constants.begin();iterisindependent){ - num_independents+= 1; - } - } - - /*Process data*/ - for(vector::iterator iter=data.begin();iterisindependent){ - num_independents+= (*iter)->M*(*iter)->N; - } - } - - /*return*/ - return num_independents; -} -/*}}}*/ -void IoModel::FillIndependents(IssmDouble* xp){/*{{{*/ - - _assert_(xp); - - /*Initialize local num ind*/ - int local_num_ind = 0; - - /*Process constants*/ - for(vector::iterator iter=constants.begin();iterisindependent){ - (*iter)->constant->GetParameterValue(&xp[local_num_ind]); - local_num_ind += 1; - } - } - - /*Process data*/ - for(vector::iterator iter=data.begin();iterisindependent){ - for(int i=0;i<(*iter)->M*(*iter)->N;i++){ - xp[local_num_ind+i] = (*iter)->data[i]; - } - local_num_ind += (*iter)->M*(*iter)->N; - } - } - - _assert_(local_num_ind == this->NumIndependents()); -} -/*}}}*/ Param* IoModel::CopyConstantObject(const char* constant_name,int param_enum){/*{{{*/ @@ -651,44 +407,4 @@ return NULL; -} -/*}}}*/ -void IoModel::StartTrace(bool trace){/*{{{*/ - - bool autodiff = false; - bool iscontrol = false; - bool keep=false; - IssmDouble gcTriggerRatio; - IssmDouble gcTriggerMaxSize; - IssmDouble obufsize; - IssmDouble lbufsize; - IssmDouble cbufsize; - IssmDouble tbufsize; - - int my_rank=IssmComm::GetRank(); - - this->FetchData(&autodiff,"md.autodiff.isautodiff"); - this->FetchData(&iscontrol,"md.inversion.iscontrol"); - - if(trace || (autodiff && !iscontrol)){ - #ifdef _HAVE_ADOLC_ - /*Retrieve parameters: */ - this->FetchData(&keep,"md.autodiff.keep"); - int keepTaylors=keep?1:0; - this->FetchData(&gcTriggerRatio,"md.autodiff.gcTriggerRatio"); - this->FetchData(&gcTriggerMaxSize,"md.autodiff.gcTriggerMaxSize"); - this->FetchData(&obufsize,"md.autodiff.obufsize"); - this->FetchData(&lbufsize,"md.autodiff.lbufsize"); - this->FetchData(&cbufsize,"md.autodiff.cbufsize"); - this->FetchData(&tbufsize,"md.autodiff.tbufsize"); - - /*Set garbage collection parameters: */ - setStoreManagerControl(reCast(gcTriggerRatio),reCast(gcTriggerMaxSize)); - - /*Start trace: */ - int skipFileDeletion=1; - trace_on(my_rank,keepTaylors,reCast(obufsize),reCast(lbufsize),reCast(cbufsize),reCast(tbufsize),skipFileDeletion); - #endif - } - } /*}}}*/ @@ -1519,4 +1235,235 @@ } /*}}}*/ +void IoModel::FetchDataToInput(Elements* elements,const char* vector_name,int input_enum,IssmDouble default_value){/*{{{*/ + + /*First, look whether it is not already loaded in this->data*/ + vector::iterator iter; + for(iter=data.begin();itername,vector_name)==0){ + _assert_(iodata->code==7); + for(int i=0;iSize();i++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + element->InputCreate(iodata->data,this,iodata->M,iodata->N,iodata->layout,input_enum,iodata->code);//we need i to index into elements. + } + return; + } + } + + /*intermediary: */ + int code,vector_layout; + IssmDouble *doublearray = NULL; + int M,N; + + /*First of, find the record for the name, and get code of data type: */ + this->SetFilePointerToData(&code, &vector_layout,vector_name); + + /*Defaulting only supported for double arrays*/ + if(code!=7) _error_(vector_name<<" is not a double array"); + + this->FetchData(&doublearray,&M,&N,vector_name); + + for(int i=0;iSize();i++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + if(!doublearray) element->AddInput(input_enum,&default_value,P0Enum); + else element->InputCreate(doublearray,this,M,N,vector_layout,input_enum,code);//we need i to index into elements. + } + + /*Free ressources*/ + xDelete(doublearray); +} +/*}}}*/ +void IoModel::FetchDataToInput(Elements* elements,const char* vector_name,int input_enum){/*{{{*/ + + /*First, look whether it is not already loaded in this->data*/ + vector::iterator iter; + for(iter=data.begin();itername,vector_name)==0){ + for(int i=0;iSize();i++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + element->InputCreate(iodata->data,this,iodata->M,iodata->N,iodata->layout,input_enum,iodata->code);//we need i to index into elements. + } + return; + } + } + + /*intermediary: */ + int i; + int code,vector_layout; + + /*variables being fetched: */ + bool boolean; + int integer; + IssmDouble scalar; + char *string = NULL; + IssmDouble *doublearray = NULL; + int M,N; + + /*First of, find the record for the name, and get code of data type: */ + this->SetFilePointerToData(&code, &vector_layout,vector_name); + + switch(code){ + case 1: //boolean constant + this->FetchData(&boolean,vector_name); + for(i=0;iSize();i++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + element->InputUpdateFromConstant(boolean,input_enum); + } + break; + case 2: //integer constant + this->FetchData(&integer,vector_name); + for(i=0;iSize();i++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + element->InputUpdateFromConstant(integer,input_enum); + } + break; + case 3: //IssmDouble constant + this->FetchData(&scalar,vector_name); + for(i=0;iSize();i++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + element->InputUpdateFromConstant(scalar,input_enum); + } + break; + case 5: //boolean vector + this->FetchData(&doublearray,&M,&N,vector_name); //we still have a doublearray, because it might include times in transient mode + if(!doublearray) _error_("\""<Size();i++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + element->InputCreate(doublearray,this,M,N,vector_layout,input_enum,code);//we need i to index into elements. + } + break; + case 6: //int vector + this->FetchData(&doublearray,&M,&N,vector_name); //we still have a doublearray, because it might include times in transient mode + if(!doublearray) _error_("\""<Size();i++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + element->InputCreate(doublearray,this,M,N,vector_layout,input_enum,code);//we need i to index into elements. + } + break; + case 7: //IssmDouble vector + this->FetchData(&doublearray,&M,&N,vector_name); + if(!doublearray) _error_("\""<Size();i++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + element->InputCreate(doublearray,this,M,N,vector_layout,input_enum,code);//we need i to index into elements. + } + break; + default: + _error_("data code " << code << " not supported yet (detected while processing \""<(doublearray); + xDelete(string); +} +/*}}}*/ +void IoModel::FetchIndependentConstant(int* pXcount,IssmPDouble* X,const char* constant_name){/*{{{*/ + + /*recover my_rank:*/ + int my_rank=IssmComm::GetRank(); + + /*recover Xcount if X is not NULL:*/ + int Xcount = 0; + if(X) Xcount=*pXcount; + + #ifdef _HAVE_ADOLC_ //cannot come here unless you are running AD mode, from DeclaredIndependents: + + /*output: */ + IssmPDouble pscalar; + IssmDouble scalar; //same as pscalar, except it's an ADOLC independent variable + int code; + + /*Set file pointer to beginning of the data: */ + fid=this->SetFilePointerToData(&code,NULL,constant_name); + if(code!=3) _error_("expecting a IssmDouble for \"" << constant_name<<"\""); + + /*We have to read a scalar from disk. First read the dimensions of the scalar, then the scalar: */ + if(my_rank==0){ + if(fread(&pscalar,sizeof(IssmPDouble),1,fid)!=1)_error_("could not read scalar "); + + /*Now, before we even broadcast this to other nodes, declare the scalar as an independent variable!. If we + *have been supplied an X vector, use it instead of what we just read: */ + if(X){ + scalar<<=X[Xcount]; + } + else{ + scalar<<=pscalar; + } + } + + ISSM_MPI_Bcast(&scalar,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + this->AddConstantIndependent(new IoConstant(scalar,constant_name)); + + /*increment offset into X vector, now that we have read 1 value:*/ + Xcount++; *pXcount=Xcount; + #endif +} +/*}}}*/ +void IoModel::FetchIndependentData(int* pXcount,IssmPDouble* X,const char* data_name){/*{{{*/ + + /*recover my_rank:*/ + int my_rank=IssmComm::GetRank(); + + /*recover Xcount if X is not NULL:*/ + int Xcount = 0; + if(X) Xcount=*pXcount; + + #ifdef _HAVE_ADOLC_ //cannot come here unless you are running AD mode, from DeclaredIndependents: + + /*Intermediaries*/ + int M,N; + IssmPDouble* buffer=NULL; //a buffer to read the data from disk + IssmDouble* matrix=NULL; //our independent variable + int code,layout; + + /*Set file pointer to beginning of the data: */ + fid=this->SetFilePointerToData(&code,&layout,data_name); + if((code!=5) && (code!=6) && (code!=7))_error_("expecting a IssmDouble, integer or boolean matrix for \"" << data_name<<"\""); + + /*We have to read a matrix from disk. First read the dimensions of the matrix, then the whole matrix: */ + /*numberofelements: */ + if(my_rank==0){ + if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows for matrix "); + } + ISSM_MPI_Bcast(&M,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + + if(my_rank==0){ + if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns for matrix "); + } + ISSM_MPI_Bcast(&N,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + + /*Now allocate matrix: */ + if(M*N){ + buffer=xNew(M*N); + matrix=xNew(M*N); + + /*Read matrix on node 0, then broadcast: */ + if(my_rank==0){ + if(fread(buffer,M*N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read matrix "); + + /*Now, before we even broadcast this to other nodes, declare the whole matrix as a independent variable! + If we have been supplied an X vector, use it instead of what we just read: */ + if(X){ + for(int i=0;i(buffer); + } + else _error_("cannot declare the independent variable \"" << data_name << "\" if it's empty!"); + + /*Add to data as independent*/ + this->AddDataIndependent(new IoData(matrix,code,layout,M,N,data_name)); + + /*increment offset into X vector, now that we have read M*N values:*/ + Xcount+=M*N; *pXcount=Xcount; + #endif +} +/*}}}*/ void IoModel::FetchMultipleData(char*** pstrings,int* pnumstrings,const char* data_name){/*{{{*/ @@ -1897,126 +1844,116 @@ } /*}}}*/ -void IoModel::FetchDataToInput(Elements* elements,const char* vector_name,int input_enum,IssmDouble default_value){/*{{{*/ - - /*First, look whether it is not already loaded in this->data*/ - vector::iterator iter; - for(iter=data.begin();itername,vector_name)==0){ - _assert_(iodata->code==7); - for(int i=0;iSize();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->InputCreate(iodata->data,this,iodata->M,iodata->N,iodata->layout,input_enum,iodata->code);//we need i to index into elements. - } +void IoModel::FillIndependents(IssmDouble* xp){/*{{{*/ + + _assert_(xp); + + /*Initialize local num ind*/ + int local_num_ind = 0; + + /*Process constants*/ + for(vector::iterator iter=constants.begin();iterisindependent){ + (*iter)->constant->GetParameterValue(&xp[local_num_ind]); + local_num_ind += 1; + } + } + + /*Process data*/ + for(vector::iterator iter=data.begin();iterisindependent){ + for(int i=0;i<(*iter)->M*(*iter)->N;i++){ + xp[local_num_ind+i] = (*iter)->data[i]; + } + local_num_ind += (*iter)->M*(*iter)->N; + } + } + + _assert_(local_num_ind == this->NumIndependents()); +} +/*}}}*/ +void IoModel::FindConstant(bool* pvalue,const char* constant_name){/*{{{*/ + + /*Intermediary*/ + vector::iterator iter; + + for(iter=constants.begin();itername,constant_name)==0){ + ioconstant->constant->GetParameterValue(pvalue); return; } } - /*intermediary: */ - int code,vector_layout; - IssmDouble *doublearray = NULL; - int M,N; - - /*First of, find the record for the name, and get code of data type: */ - this->SetFilePointerToData(&code, &vector_layout,vector_name); - - /*Defaulting only supported for double arrays*/ - if(code!=7) _error_(vector_name<<" is not a double array"); - - this->FetchData(&doublearray,&M,&N,vector_name); - - for(int i=0;iSize();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - if(!doublearray) element->AddInput(input_enum,&default_value,P0Enum); - else element->InputCreate(doublearray,this,M,N,vector_layout,input_enum,code);//we need i to index into elements. - } - - /*Free ressources*/ - xDelete(doublearray); -} -/*}}}*/ -void IoModel::FetchDataToInput(Elements* elements,const char* vector_name,int input_enum){/*{{{*/ - - /*First, look whether it is not already loaded in this->data*/ - vector::iterator iter; - for(iter=data.begin();itername,vector_name)==0){ - for(int i=0;iSize();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->InputCreate(iodata->data,this,iodata->M,iodata->N,iodata->layout,input_enum,iodata->code);//we need i to index into elements. - } + for(vector::iterator iter=constants.begin();iterconstant->Echo(); + _error_("Could not find constant \""<::iterator iter; + + for(iter=constants.begin();itername,constant_name)==0){ + ioconstant->constant->GetParameterValue(pvalue); return; } } - /*intermediary: */ - int i; - int code,vector_layout; - - /*variables being fetched: */ - bool boolean; - int integer; - IssmDouble scalar; - char *string = NULL; - IssmDouble *doublearray = NULL; - int M,N; - - /*First of, find the record for the name, and get code of data type: */ - this->SetFilePointerToData(&code, &vector_layout,vector_name); - - switch(code){ - case 1: //boolean constant - this->FetchData(&boolean,vector_name); - for(i=0;iSize();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->InputUpdateFromConstant(boolean,input_enum); - } - break; - case 2: //integer constant - this->FetchData(&integer,vector_name); - for(i=0;iSize();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->InputUpdateFromConstant(integer,input_enum); - } - break; - case 3: //IssmDouble constant - this->FetchData(&scalar,vector_name); - for(i=0;iSize();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->InputUpdateFromConstant(scalar,input_enum); - } - break; - case 5: //boolean vector - this->FetchData(&doublearray,&M,&N,vector_name); //we still have a doublearray, because it might include times in transient mode - if(!doublearray) _error_("\""<Size();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->InputCreate(doublearray,this,M,N,vector_layout,input_enum,code);//we need i to index into elements. - } - break; - case 6: //int vector - this->FetchData(&doublearray,&M,&N,vector_name); //we still have a doublearray, because it might include times in transient mode - if(!doublearray) _error_("\""<Size();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->InputCreate(doublearray,this,M,N,vector_layout,input_enum,code);//we need i to index into elements. - } - break; - case 7: //IssmDouble vector - this->FetchData(&doublearray,&M,&N,vector_name); - if(!doublearray) _error_("\""<Size();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->InputCreate(doublearray,this,M,N,vector_layout,input_enum,code);//we need i to index into elements. - } - break; - default: - _error_("data code " << code << " not supported yet (detected while processing \""<(doublearray); - xDelete(string); + _error_("Could not find constant \""<::iterator iter; + + for(iter=constants.begin();itername,constant_name)==0){ + ioconstant->constant->GetParameterValue(pvalue); + return; + } + } + + _error_("Could not find constant \""<::iterator iter; + + for(iter=constants.begin();itername,constant_name)==0){ + ioconstant->constant->GetParameterValue(pvalue); + return; + } + } + + _error_("Could not find constant \""<::iterator iter; + + for(iter=constants.begin();itername,constant_name)==0){ + ioconstant->constant->GetParameterValue(pvalue,psize); + return; + } + } + + _error_("Could not find constant \""<(record_name); - break; - } - if(record_name_size<3 || record_name_size>80){ - _error_("error while looking in binary file. Found a string of size "<(record_name_size+1); - record_name[record_name_size]='\0'; - - /*Read record_name: */ - if(fread(record_name,record_name_size*sizeof(char),1,fid)==0){ - /*we have reached the end of the file. break: */ - found=0; - xDelete(record_name); - break; - } - if(strncmp(record_name,mddot,3)!=0){ - _error_("error while reading binary file: record does not start with \"md.\": "<(record_name); - break; - } - else{ - /*This is not the correct string, read the record length, and use it to skip this record: */ - if(fread(&record_length,sizeof(int),1,fid)!=1) _error_("Could not read record_length"); - /*skip: */ - fseek(fid,record_length,SEEK_CUR); - xDelete(record_name); - } - } - } - ISSM_MPI_Bcast(&found,1,ISSM_MPI_INT,0,IssmComm::GetComm()); - if(!found)_error_("could not find data with name \"" << data_name << "\" in binary file"); - - /*Broadcast code and vector type: */ - ISSM_MPI_Bcast(&record_code,1,ISSM_MPI_INT,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(&vector_type,1,ISSM_MPI_INT,0,IssmComm::GetComm()); - - /*Assign output pointers:*/ - *pcode=record_code; - if(pvector_type)*pvector_type=vector_type; - - return fid; +int IoModel::NumIndependents(void){/*{{{*/ + + /*Initialize output*/ + int num_independents = 0; + + /*Process constants*/ + for(vector::iterator iter=constants.begin();iterisindependent){ + num_independents+= 1; + } + } + + /*Process data*/ + for(vector::iterator iter=data.begin();iterisindependent){ + num_independents+= (*iter)->M*(*iter)->N; + } + } + + /*return*/ + return num_independents; } /*}}}*/ @@ -2284,2 +2158,128 @@ } /*}}}*/ +FILE* IoModel::SetFilePointerToData(int* pcode,int* pvector_type,const char* data_name){/*{{{*/ + + int my_rank; + + int found = 0; + const char* mddot = "md."; + char* record_name = NULL; + int record_name_size; + int record_length; + int record_code; //1 to 7 number + int vector_type = 0; //nodal or elementary + + /*recover my_rank:*/ + my_rank=IssmComm::GetRank(); + _assert_(strncmp(data_name,mddot,3)==0); + + /*Go find in the binary file, the position of the data we want to fetch: */ + if(my_rank==0){ + + /*First set FILE* position to the beginning of the file: */ + fseek(fid,0,SEEK_SET); + + /*Now march through file looking for the correct data identifier: */ + for(;;){ + /*Read size of first string name: */ + if(fread(&record_name_size,sizeof(int),1,fid)==0){ + /*we have reached the end of the file. break: */ + xDelete(record_name); + break; + } + if(record_name_size<3 || record_name_size>80){ + _error_("error while looking in binary file. Found a string of size "<(record_name_size+1); + record_name[record_name_size]='\0'; + + /*Read record_name: */ + if(fread(record_name,record_name_size*sizeof(char),1,fid)==0){ + /*we have reached the end of the file. break: */ + found=0; + xDelete(record_name); + break; + } + if(strncmp(record_name,mddot,3)!=0){ + _error_("error while reading binary file: record does not start with \"md.\": "<(record_name); + break; + } + else{ + /*This is not the correct string, read the record length, and use it to skip this record: */ + if(fread(&record_length,sizeof(int),1,fid)!=1) _error_("Could not read record_length"); + /*skip: */ + fseek(fid,record_length,SEEK_CUR); + xDelete(record_name); + } + } + } + ISSM_MPI_Bcast(&found,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + if(!found)_error_("could not find data with name \"" << data_name << "\" in binary file"); + + /*Broadcast code and vector type: */ + ISSM_MPI_Bcast(&record_code,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&vector_type,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + + /*Assign output pointers:*/ + *pcode=record_code; + if(pvector_type)*pvector_type=vector_type; + + return fid; +} +/*}}}*/ +void IoModel::StartTrace(bool trace){/*{{{*/ + + bool autodiff = false; + bool iscontrol = false; + bool keep=false; + IssmDouble gcTriggerRatio; + IssmDouble gcTriggerMaxSize; + IssmDouble obufsize; + IssmDouble lbufsize; + IssmDouble cbufsize; + IssmDouble tbufsize; + + int my_rank=IssmComm::GetRank(); + + this->FetchData(&autodiff,"md.autodiff.isautodiff"); + this->FetchData(&iscontrol,"md.inversion.iscontrol"); + + if(trace || (autodiff && !iscontrol)){ + #ifdef _HAVE_ADOLC_ + /*Retrieve parameters: */ + this->FetchData(&keep,"md.autodiff.keep"); + int keepTaylors=keep?1:0; + this->FetchData(&gcTriggerRatio,"md.autodiff.gcTriggerRatio"); + this->FetchData(&gcTriggerMaxSize,"md.autodiff.gcTriggerMaxSize"); + this->FetchData(&obufsize,"md.autodiff.obufsize"); + this->FetchData(&lbufsize,"md.autodiff.lbufsize"); + this->FetchData(&cbufsize,"md.autodiff.cbufsize"); + this->FetchData(&tbufsize,"md.autodiff.tbufsize"); + + /*Set garbage collection parameters: */ + setStoreManagerControl(reCast(gcTriggerRatio),reCast(gcTriggerMaxSize)); + + /*Start trace: */ + int skipFileDeletion=1; + trace_on(my_rank,keepTaylors,reCast(obufsize),reCast(lbufsize),reCast(cbufsize),reCast(tbufsize),skipFileDeletion); + #endif + } + +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/IoModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/IoModel.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/IoModel.h (revision 20810) @@ -18,7 +18,7 @@ class IoConstant { /*holds single IssmDouble, int, bool and char from input*/ public: + Param* constant; bool isindependent; char* name; - Param* constant; ~IoConstant(); @@ -33,10 +33,10 @@ class IoData { /*holds temporary data (array), memory intensive*/ public: + int code; + IssmDouble* data; bool isindependent; + int layout; + int M,N; char* name; - int M,N; - int code; - int layout; - IssmDouble* data; ~IoData(); @@ -67,24 +67,24 @@ /*Mesh properties and connectivity tables*/ + int domaindim; int domaintype; - int domaindim; - int meshelementtype; - int numberofvertices; - int numberofelements; - int numberoffaces; - int numberofedges; - int facescols; int *elements; - int *faces; int *edges; int *elementtoedgeconnectivity; int *elementtofaceconnectivity; + int *faces; + int facescols; + int meshelementtype; + int *numbernodetoelementconnectivity; + int numberofedges; + int numberofelements; + int numberoffaces; + int numberofvertices; int *singlenodetoelementconnectivity; - int *numbernodetoelementconnectivity; /*Data to synchronize through low level object drivers: */ + int constraintcounter; //keep track of how many constraints are being created in each analysis + int loadcounter; //keep track of how many loads are being created in each analysis int nodecounter; //keep track of how many nodes are being created in each analysis - int loadcounter; //keep track of how many loads are being created in each analysis - int constraintcounter; //keep track of how many constraints are being created in each analysis /*Methods*/ @@ -98,4 +98,7 @@ void AddData(IoData* data_in); void AddDataIndependent(IoData* data_in); + void FetchIndependentConstant(int* pXcount,IssmPDouble* X,const char* name); + void FetchIndependentData(int* pXcount,IssmPDouble* X,const char* name); + void FillIndependents(IssmDouble* xp); void FindConstant(bool* pvalue,const char* constant_name); void FindConstant(int* pvalue,const char* constant_name); @@ -103,8 +106,5 @@ void FindConstant(char **pvalue,const char* constant_name); void FindConstant(char ***pvalue,int* psize,const char* constant_name); - void FetchIndependentConstant(int* pXcount,IssmPDouble* X,const char* name); - void FetchIndependentData(int* pXcount,IssmPDouble* X,const char* name); int NumIndependents(); - void FillIndependents(IssmDouble* xp); /*Input/Output*/ @@ -112,4 +112,5 @@ Param *CopyConstantObject(const char* constant_name,int param_enum); IssmDouble *Data(const char* data_name); + void DeclareIndependents(bool trace,IssmPDouble* X); void DeleteData(int num,...); void DeleteData(IssmDouble* vector,const char* data_name); @@ -123,4 +124,9 @@ void FetchData(IssmDouble** pscalarmatrix,int* pM,int* pN,const char* data_name); void FetchData(IssmDouble*** pmatrixarray,int** pmdims,int** pndims, int* pnumrecords,const char* data_name); + void FetchData(Option **poption,int index); + void FetchData(int num,...); + void FetchDataToInput(Elements* elements,const char* vector_name,int input_enum); + void FetchDataToInput(Elements* elements,const char* vector_name,int input_enum,IssmDouble default_value); + void FetchIndependent(const char* dependent_name); void FetchMultipleData(char*** pstringarray,int* pnumstrings,const char* data_name); void FetchMultipleData(IssmDouble*** pmatrixarray,int** pmdims,int** pndims, int* pnumrecords,const char* data_name); @@ -128,14 +134,8 @@ void FetchMultipleData(int** pvector, int* pnum_instances,const char* data_name); void FetchMultipleData(IssmDouble** pvector, int* pnum_instances,const char* data_name); - void FetchData(Option **poption,int index); - void FetchData(int num,...); - void FetchDataToInput(Elements* elements,const char* vector_name,int input_enum); - void FetchDataToInput(Elements* elements,const char* vector_name,int input_enum,IssmDouble default_value); void LastIndex(int *pindex); + fpos_t* SetFilePointersToData(int** pcodes,int** pvector_types, int* pnum_instances, const char* data_name); FILE* SetFilePointerToData(int* pcode,int* pvector_type, const char* data_name); - fpos_t* SetFilePointersToData(int** pcodes,int** pvector_types, int* pnum_instances, const char* data_name); - void DeclareIndependents(bool trace,IssmPDouble* X); void StartTrace(bool trace); - void FetchIndependent(const char* dependent_name); }; Index: /issm/trunk-jpl/src/c/classes/Masscon.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Masscon.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Masscon.h (revision 20810) @@ -54,4 +54,13 @@ /*}}}*/ /*Object virtual function resolutoin: */ + Object* copy() {/*{{{*/ + Masscon* mf = new Masscon(this->name,this->definitionenum,this->levelset,this->M); + return (Object*) mf; + } + /*}}}*/ + void DeepEcho(void){/*{{{*/ + this->Echo(); + } + /*}}}*/ void Echo(void){/*{{{*/ _printf_(" Masscon: " << this->name << " " << this->definitionenum << "\n"); @@ -60,19 +69,6 @@ } /*}}}*/ - void DeepEcho(void){/*{{{*/ - this->Echo(); - } - /*}}}*/ int Id(void){/*{{{*/ return -1; - } - /*}}}*/ - int ObjectEnum(void){/*{{{*/ - return MassconEnum; - } - /*}}}*/ - Object* copy() {/*{{{*/ - Masscon* mf = new Masscon(this->name,this->definitionenum,this->levelset,this->M); - return (Object*) mf; } /*}}}*/ @@ -81,5 +77,14 @@ } /*}}}*/ + int ObjectEnum(void){/*{{{*/ + return MassconEnum; + } + /*}}}*/ /*Definition virtual function resolutoin: */ + int DefinitionEnum(){/*{{{*/ + + return this->definitionenum; + } + /*}}}*/ char* Name(){/*{{{*/ @@ -88,9 +93,4 @@ return name2; - } - /*}}}*/ - int DefinitionEnum(){/*{{{*/ - - return this->definitionenum; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Massconaxpby.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Massconaxpby.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Massconaxpby.h (revision 20810) @@ -63,4 +63,13 @@ /*}}}*/ /*Object virtual function resolutoin: */ + Object* copy() {/*{{{*/ + Massconaxpby* mf = new Massconaxpby(this->name,this->definitionenum,this->namex,this->namey, this->alpha, this->beta); + return (Object*) mf; + } + /*}}}*/ + void DeepEcho(void){/*{{{*/ + this->Echo(); + } + /*}}}*/ void Echo(void){/*{{{*/ _printf_(" Massconaxpby: " << this->name << " " << this->definitionenum << "\n"); @@ -71,19 +80,6 @@ } /*}}}*/ - void DeepEcho(void){/*{{{*/ - this->Echo(); - } - /*}}}*/ int Id(void){/*{{{*/ return -1; - } - /*}}}*/ - int ObjectEnum(void){/*{{{*/ - return MassconaxpbyEnum; - } - /*}}}*/ - Object* copy() {/*{{{*/ - Massconaxpby* mf = new Massconaxpby(this->name,this->definitionenum,this->namex,this->namey, this->alpha, this->beta); - return (Object*) mf; } /*}}}*/ @@ -92,5 +88,14 @@ } /*}}}*/ + int ObjectEnum(void){/*{{{*/ + return MassconaxpbyEnum; + } + /*}}}*/ /*Definition virtual function resolutoin: */ + int DefinitionEnum(){/*{{{*/ + + return this->definitionenum; + } + /*}}}*/ char* Name(){/*{{{*/ @@ -99,9 +104,4 @@ return name2; - } - /*}}}*/ - int DefinitionEnum(){/*{{{*/ - - return this->definitionenum; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Massfluxatgate.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Massfluxatgate.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Massfluxatgate.h (revision 20810) @@ -94,4 +94,12 @@ /*Object virtual function resolutoin: */ + Object* copy() {/*{{{*/ + return new Massfluxatgate(this->name,this->definitionenum,this->numsegments,this->x1,this->y1,this->x2,this->y2,this->elements); + } + /*}}}*/ + void DeepEcho(void){/*{{{*/ + this->Echo(); + } + /*}}}*/ void Echo(void){/*{{{*/ _printf_(" Massfluxatgate: " << name << " " << this->definitionenum << "\n"); @@ -105,8 +113,4 @@ } /*}}}*/ - void DeepEcho(void){/*{{{*/ - this->Echo(); - } - /*}}}*/ int Id(void){/*{{{*/ return -1; @@ -117,8 +121,4 @@ } /*}}}*/ - Object* copy() {/*{{{*/ - return new Massfluxatgate(this->name,this->definitionenum,this->numsegments,this->x1,this->y1,this->x2,this->y2,this->elements); - } - /*}}}*/ void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){/*{{{*/ _error_("not implemented yet!"); @@ -126,4 +126,9 @@ /*}}}*/ /*Definition virtual function resolutoin: */ + int DefinitionEnum(){/*{{{*/ + + return this->definitionenum; + } + /*}}}*/ char* Name(){/*{{{*/ @@ -132,9 +137,4 @@ return name2; - } - /*}}}*/ - int DefinitionEnum(){/*{{{*/ - - return this->definitionenum; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Misfit.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Misfit.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Misfit.h (revision 20810) @@ -25,13 +25,13 @@ int definitionenum; + bool local; + int model_enum; char* name; - int model_enum; int observation_enum; + char* timeinterpolation; int weights_enum; - char* timeinterpolation; - bool local; + int lock; // if lock is on, we just return the value stored in "misfit". this is used so we don't compute misfit past the final_time IssmDouble misfit; //value carried over in time. - int lock; // if lock is on, we just return the value stored in "misfit". this is used so we don't compute misfit past the final_time /*Misfit constructors, destructors :*/ @@ -76,4 +76,15 @@ /*}}}*/ /*Object virtual function resolutoin: */ + Object* copy() {/*{{{*/ + Misfit* mf = new Misfit(this->name,this->definitionenum, this->model_enum,this->observation_enum,this->timeinterpolation,this->local,this->weights_enum); + mf->misfit=this->misfit; + mf->lock=this->lock; + return (Object*) mf; + } + /*}}}*/ + void DeepEcho(void){/*{{{*/ + this->Echo(); + } + /*}}}*/ void Echo(void){/*{{{*/ _printf_(" Misfit: " << name << " " << this->definitionenum << "\n"); @@ -85,21 +96,6 @@ } /*}}}*/ - void DeepEcho(void){/*{{{*/ - this->Echo(); - } - /*}}}*/ int Id(void){/*{{{*/ return -1; - } - /*}}}*/ - int ObjectEnum(void){/*{{{*/ - return MisfitEnum; - } - /*}}}*/ - Object* copy() {/*{{{*/ - Misfit* mf = new Misfit(this->name,this->definitionenum, this->model_enum,this->observation_enum,this->timeinterpolation,this->local,this->weights_enum); - mf->misfit=this->misfit; - mf->lock=this->lock; - return (Object*) mf; } /*}}}*/ @@ -108,5 +104,14 @@ } /*}}}*/ + int ObjectEnum(void){/*{{{*/ + return MisfitEnum; + } + /*}}}*/ /*Definition virtual function resolutoin: */ + int DefinitionEnum(){/*{{{*/ + + return this->definitionenum; + } + /*}}}*/ char* Name(){/*{{{*/ @@ -115,9 +120,4 @@ return name2; - } - /*}}}*/ - int DefinitionEnum(){/*{{{*/ - - return this->definitionenum; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Nodalvalue.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Nodalvalue.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Nodalvalue.h (revision 20810) @@ -27,6 +27,6 @@ int definitionenum; + int model_enum; char* name; - int model_enum; int node; @@ -56,4 +56,13 @@ /*}}}*/ /*Object virtual function resolutoin: */ + Object* copy() {/*{{{*/ + Nodalvalue* mf = new Nodalvalue(this->name,this->definitionenum, this->model_enum,this->node); + return (Object*) mf; + } + /*}}}*/ + void DeepEcho(void){/*{{{*/ + this->Echo(); + } + /*}}}*/ void Echo(void){/*{{{*/ _printf_(" Nodalvalue: " << name << " " << this->definitionenum << "\n"); @@ -62,19 +71,6 @@ } /*}}}*/ - void DeepEcho(void){/*{{{*/ - this->Echo(); - } - /*}}}*/ int Id(void){/*{{{*/ return -1; - } - /*}}}*/ - int ObjectEnum(void){/*{{{*/ - return NodalvalueEnum; - } - /*}}}*/ - Object* copy() {/*{{{*/ - Nodalvalue* mf = new Nodalvalue(this->name,this->definitionenum, this->model_enum,this->node); - return (Object*) mf; } /*}}}*/ @@ -83,5 +79,14 @@ } /*}}}*/ + int ObjectEnum(void){/*{{{*/ + return NodalvalueEnum; + } + /*}}}*/ /*Definition virtual function resolutoin: */ + int DefinitionEnum(){/*{{{*/ + + return this->definitionenum; + } + /*}}}*/ char* Name(){/*{{{*/ @@ -90,9 +95,4 @@ return name2; - } - /*}}}*/ - int DefinitionEnum(){/*{{{*/ - - return this->definitionenum; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Node.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Node.cpp (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Node.cpp (revision 20810) @@ -165,4 +165,16 @@ /*Object virtual functions definitions:*/ +void Node::DeepEcho(void){/*{{{*/ + + _printf_("Node:\n"); + _printf_(" id: " << id << "\n"); + _printf_(" sid: " << sid << "\n"); + _printf_(" analysis_enum: " << EnumToStringx(analysis_enum) << "\n"); + _printf_(" approximation: " << EnumToStringx(approximation) << "\n"); + _printf_(" indexingupdate: " << indexingupdate << "\n"); + indexing.DeepEcho(); + +} +/*}}}*/ void Node::Echo(void){/*{{{*/ @@ -177,16 +189,4 @@ } /*}}}*/ -void Node::DeepEcho(void){/*{{{*/ - - _printf_("Node:\n"); - _printf_(" id: " << id << "\n"); - _printf_(" sid: " << sid << "\n"); - _printf_(" analysis_enum: " << EnumToStringx(analysis_enum) << "\n"); - _printf_(" approximation: " << EnumToStringx(approximation) << "\n"); - _printf_(" indexingupdate: " << indexingupdate << "\n"); - indexing.DeepEcho(); - -} -/*}}}*/ int Node::Id(void){ return id; }/*{{{*/ /*}}}*/ @@ -199,4 +199,11 @@ /*Node management:*/ +void Node::GetCoordinateSystem(IssmDouble* coord_system_out){/*{{{*/ + + /*Copy coord_system*/ + for(int k=0;k<3;k++) for(int l=0;l<3;l++) coord_system_out[3*k+l]=this->coord_system[k][l]; + +} +/*}}}*/ int Node::GetDof(int dofindex,int setenum){/*{{{*/ @@ -380,19 +387,4 @@ } /*}}}*/ -int Node::Sid(void){/*{{{*/ - return sid; -} -/*}}}*/ -int Node::Lid(void){/*{{{*/ - return lid; -} -/*}}}*/ -void Node::GetCoordinateSystem(IssmDouble* coord_system_out){/*{{{*/ - - /*Copy coord_system*/ - for(int k=0;k<3;k++) for(int l=0;l<3;l++) coord_system_out[3*k+l]=this->coord_system[k][l]; - -} -/*}}}*/ bool Node::InAnalysis(int in_analysis_enum){/*{{{*/ if (in_analysis_enum==this->analysis_enum) return true; @@ -400,6 +392,23 @@ } /*}}}*/ +int Node::Lid(void){/*{{{*/ + return lid; +} +/*}}}*/ +int Node::Sid(void){/*{{{*/ + return sid; +} +/*}}}*/ /*Node numerics:*/ +void Node::Activate(void){/*{{{*/ + + if(!IsActive() && !this->indexing.freeze){ + this->indexingupdate = true; + indexing.Activate(); + } + +} +/*}}}*/ void Node::ApplyConstraint(int dof,IssmDouble value){/*{{{*/ @@ -408,23 +417,4 @@ DofInSSet(dof); this->indexing.svalues[dof]=value; -} -/*}}}*/ -bool Node::RequiresDofReindexing(void){/*{{{*/ - - return this->indexingupdate; - -} -/*}}}*/ -void Node::ReindexingDone(void){/*{{{*/ - - this->indexingupdate = false; - -} -/*}}}*/ -void Node::RelaxConstraint(int dof){/*{{{*/ - - /*Dof should be added to the f-set, and taken out of the s-set:*/ - DofInFSet(dof); - this->indexing.svalues[dof]=0.; } /*}}}*/ @@ -456,4 +446,27 @@ } /*}}}*/ +void Node::Deactivate(void){/*{{{*/ + + if(IsActive() && !this->indexing.freeze){ + this->indexingupdate = true; + indexing.Deactivate(); + } + +} +/*}}}*/ +void Node::DofInFSet(int dof){/*{{{*/ + + /*Put dof for this node into the f set (ie, this dof will NOT be constrained + * to a fixed value during computations. */ + _assert_(dofindexing.gsize); + + if(this->indexing.f_set[dof] == 0){ + if(this->indexing.freeze) _error_("Cannot change dof of frozen node"); + this->indexingupdate = true; + this->indexing.f_set[dof]=1; + this->indexing.s_set[dof]=0; + } +} +/*}}}*/ void Node::DofInSSet(int dof){/*{{{*/ @@ -470,18 +483,4 @@ } /*}}}*/ -void Node::DofInFSet(int dof){/*{{{*/ - - /*Put dof for this node into the f set (ie, this dof will NOT be constrained - * to a fixed value during computations. */ - _assert_(dofindexing.gsize); - - if(this->indexing.f_set[dof] == 0){ - if(this->indexing.freeze) _error_("Cannot change dof of frozen node"); - this->indexingupdate = true; - this->indexing.f_set[dof]=1; - this->indexing.s_set[dof]=0; - } -} -/*}}}*/ void Node::FreezeDof(int dof){/*{{{*/ @@ -489,35 +488,4 @@ //FIXME: for now we don't want this element to change so we use freeze this->indexing.freeze =true; - -} -/*}}}*/ -void Node::Deactivate(void){/*{{{*/ - - if(IsActive() && !this->indexing.freeze){ - this->indexingupdate = true; - indexing.Deactivate(); - } - -} -/*}}}*/ -void Node::HardDeactivate(void){/*{{{*/ - - this->indexing.Deactivate(); - this->indexing.freeze =true; - -} -/*}}}*/ -void Node::Activate(void){/*{{{*/ - - if(!IsActive() && !this->indexing.freeze){ - this->indexingupdate = true; - indexing.Activate(); - } - -} -/*}}}*/ -bool Node::IsActive(void){/*{{{*/ - - return indexing.active; } @@ -575,7 +543,39 @@ } /*}}}*/ +void Node::HardDeactivate(void){/*{{{*/ + + this->indexing.Deactivate(); + this->indexing.freeze =true; + +} +/*}}}*/ +bool Node::IsActive(void){/*{{{*/ + + return indexing.active; + +} +/*}}}*/ int Node::IsClone(){/*{{{*/ return indexing.clone; + +} +/*}}}*/ +void Node::ReindexingDone(void){/*{{{*/ + + this->indexingupdate = false; + +} +/*}}}*/ +void Node::RelaxConstraint(int dof){/*{{{*/ + + /*Dof should be added to the f-set, and taken out of the s-set:*/ + DofInFSet(dof); + this->indexing.svalues[dof]=0.; +} +/*}}}*/ +bool Node::RequiresDofReindexing(void){/*{{{*/ + + return this->indexingupdate; } @@ -738,4 +738,21 @@ } else _error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!"); +} +/*}}}*/ +void Node::SetClone(int* minranks){/*{{{*/ + + int my_rank; + + /*recover my_rank:*/ + my_rank=IssmComm::GetRank(); + + if (minranks[sid]==my_rank){ + indexing.clone=false; + } + else{ + /*!there is a cpu with lower rank that has the same node, + therefore, I am a clone*/ + indexing.clone=true; + } } /*}}}*/ @@ -788,23 +805,43 @@ } /*}}}*/ -void Node::SetClone(int* minranks){/*{{{*/ - - int my_rank; - - /*recover my_rank:*/ - my_rank=IssmComm::GetRank(); - - if (minranks[sid]==my_rank){ - indexing.clone=false; - } - else{ - /*!there is a cpu with lower rank that has the same node, - therefore, I am a clone*/ - indexing.clone=true; - } -} -/*}}}*/ /*Methods inherent to Node: */ +int* GetGlobalDofList(Node** nodes,int numnodes,int setenum,int approximation){/*{{{*/ + + int i,numdof,count; + int* ndof_list=NULL; + int *doflist = NULL; + + if(numnodes){ + + /*Allocate:*/ + ndof_list=xNew(numnodes); + + /*First, figure out size of doflist: */ + numdof=0; + for(i=0;iGetNumberOfDofs(approximation,setenum); + numdof+=ndof_list[i]; + } + + if(numdof){ + /*Allocate: */ + doflist=xNew(numdof); + + /*Populate: */ + count=0; + for(i=0;iGetDofList(&doflist[count],approximation,setenum); + count+=ndof_list[i]; + } + } + else doflist=NULL; + } + /*Free ressources:*/ + xDelete(ndof_list); + + return doflist; +} +/*}}}*/ int* GetLocalDofList(Node** nodes,int numnodes,int setenum,int approximation){ /*{{{*/ @@ -864,41 +901,4 @@ } /*}}}*/ -int* GetGlobalDofList(Node** nodes,int numnodes,int setenum,int approximation){/*{{{*/ - - int i,numdof,count; - int* ndof_list=NULL; - int *doflist = NULL; - - if(numnodes){ - - /*Allocate:*/ - ndof_list=xNew(numnodes); - - /*First, figure out size of doflist: */ - numdof=0; - for(i=0;iGetNumberOfDofs(approximation,setenum); - numdof+=ndof_list[i]; - } - - if(numdof){ - /*Allocate: */ - doflist=xNew(numdof); - - /*Populate: */ - count=0; - for(i=0;iGetDofList(&doflist[count],approximation,setenum); - count+=ndof_list[i]; - } - } - else doflist=NULL; - } - /*Free ressources:*/ - xDelete(ndof_list); - - return doflist; -} -/*}}}*/ int GetNumberOfDofs(Node** nodes,int numnodes,int setenum,int approximation){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Node.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Node.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Node.h (revision 20810) @@ -34,8 +34,8 @@ int lid; // "local" id (rank of this node in current partition) + int analysis_enum; + IssmDouble coord_system[3][3]; bool indexingupdate; DofIndexing indexing; - int analysis_enum; - IssmDouble coord_system[3][3]; /*Node constructors, destructors*/ @@ -45,49 +45,49 @@ /*Object virtual functions definitions:*/ + Object *copy(); + void DeepEcho(); void Echo(); - void DeepEcho(); int Id(); + void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction); int ObjectEnum(); - Object *copy(); - void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction); /*Node numerical routines*/ + void Activate(void); + void ApplyConstraint(int dof,IssmDouble value); void CreateNodalConstraints(Vector* ys); - void SetCurrentConfiguration(DataSet* nodes,Vertices* vertices); - int Sid(void); - int Lid(void); + void Deactivate(void); + void DistributeDofs(int* pdofcount,int setenum); + void DofInFSet(int dof); + void DofInSSet(int dof); + void FreezeDof(int dof); + int GetApproximation(); void GetCoordinateSystem(IssmDouble* coord_system_out); - bool InAnalysis(int analysis_enum); - int GetApproximation(); - int GetNumberOfDofs(int approximation_enum,int setenum); - int IsClone(); - void ApplyConstraint(int dof,IssmDouble value); - void RelaxConstraint(int dof); - void DofInSSet(int dof); - void DofInFSet(int dof); int GetDof(int dofindex,int setenum); void GetDofList(int* poutdoflist,int approximation_enum,int setenum); void GetLocalDofList(int* poutdoflist,int approximation_enum,int setenum); - void FreezeDof(int dof); + int GetNumberOfDofs(int approximation_enum,int setenum); + void HardDeactivate(void); + bool InAnalysis(int analysis_enum); bool IsActive(void); - void Activate(void); - void Deactivate(void); - void HardDeactivate(void); - void ReindexingDone(void); - bool RequiresDofReindexing(void); + int IsClone(); int IsFloating(); int IsGrounded(); + int Lid(void); + void OffsetDofs(int dofcount,int setenum); + void ReindexingDone(void); + void RelaxConstraint(int dof); + bool RequiresDofReindexing(void); + void SetClone(int* minranks); + void SetCurrentConfiguration(DataSet* nodes,Vertices* vertices); + void ShowTrueDofs(int* truerows,int ncols,int setenum); + int Sid(void); + void UpdateCloneDofs(int* alltruerows,int ncols,int setenum); void VecMerge(Vector* ug, IssmDouble* vector_serial,int setenum); void VecReduce(Vector* vector, IssmDouble* ug_serial,int setnum); - void DistributeDofs(int* pdofcount,int setenum); - void OffsetDofs(int dofcount,int setenum); - void ShowTrueDofs(int* truerows,int ncols,int setenum); - void UpdateCloneDofs(int* alltruerows,int ncols,int setenum); - void SetClone(int* minranks); }; /*Methods inherent to Node: */ +int* GetGlobalDofList(Node** nodes,int numnodes,int setenum,int approximation); int* GetLocalDofList(Node** nodes,int numnodes,int setenum,int approximation); -int* GetGlobalDofList(Node** nodes,int numnodes,int setenum,int approximation); int GetNumberOfDofs(Node** nodes,int numnodes,int setenum,int approximation); Index: /issm/trunk-jpl/src/c/classes/Nodes.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Nodes.cpp (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Nodes.cpp (revision 20810) @@ -164,29 +164,4 @@ } /*}}}*/ -int Nodes::MaxNumDofs(int analysis_type,int setenum){/*{{{*/ - - int max=0; - int allmax,numdofs; - - /*Now go through all nodes, and get how many dofs they own, unless they are clone nodes: */ - for(int i=0;iSize();i++){ - - Node* node=xDynamicCast(this->GetObjectByOffset(i)); - - /*Check that this node corresponds to our analysis currently being carried out: */ - if (node->InAnalysis(analysis_type)){ - - numdofs=node->GetNumberOfDofs(NoneApproximationEnum,setenum); - if(numdofs>max)max=numdofs; - } - } - - /*Grab max of all cpus: */ - ISSM_MPI_Allreduce((void*)&max,(void*)&allmax,1,ISSM_MPI_INT,ISSM_MPI_MAX,IssmComm::GetComm()); - max=allmax; - - return max; -} -/*}}}*/ int Nodes::MaximumId(){/*{{{*/ @@ -219,4 +194,29 @@ } /*}}}*/ +int Nodes::MaxNumDofs(int analysis_type,int setenum){/*{{{*/ + + int max=0; + int allmax,numdofs; + + /*Now go through all nodes, and get how many dofs they own, unless they are clone nodes: */ + for(int i=0;iSize();i++){ + + Node* node=xDynamicCast(this->GetObjectByOffset(i)); + + /*Check that this node corresponds to our analysis currently being carried out: */ + if (node->InAnalysis(analysis_type)){ + + numdofs=node->GetNumberOfDofs(NoneApproximationEnum,setenum); + if(numdofs>max)max=numdofs; + } + } + + /*Grab max of all cpus: */ + ISSM_MPI_Allreduce((void*)&max,(void*)&allmax,1,ISSM_MPI_INT,ISSM_MPI_MAX,IssmComm::GetComm()); + max=allmax; + + return max; +} +/*}}}*/ int Nodes::NumberOfDofs(int analysis_type,int setenum){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Nodes.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Nodes.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Nodes.h (revision 20810) @@ -27,7 +27,6 @@ void DistributeDofs(int analysis_type,int SETENUM); void FlagClones(int analysis_type); - bool RequiresDofReindexing(int analysis_type); + int MaximumId(void); int MaxNumDofs(int analysis_type,int setenum); - int MaximumId(void); int NumberOfDofs(int analysis_type,int setenum); int NumberOfDofsLocal(int analysis_type,int setenum); @@ -35,4 +34,5 @@ int NumberOfNodes(void); void Ranks(int* ranks,int analysis_type); + bool RequiresDofReindexing(int analysis_type); }; Index: /issm/trunk-jpl/src/c/classes/Profiler.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Profiler.cpp (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Profiler.cpp (revision 20810) @@ -41,4 +41,10 @@ /*Object virtual functions definitions:*/ +void Profiler::DeepEcho(void){/*{{{*/ + + this->Echo(); + +} +/*}}}*/ void Profiler::Echo(void){/*{{{*/ @@ -49,17 +55,5 @@ } /*}}}*/ -void Profiler::DeepEcho(void){/*{{{*/ - - this->Echo(); - -} -/*}}}*/ int Profiler::Id(void){ return -1; }/*{{{*/ -/*}}}*/ -int Profiler::ObjectEnum(void){/*{{{*/ - - return ProfilerEnum; - -} /*}}}*/ void Profiler::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/ @@ -76,6 +70,78 @@ } /*}}}*/ +int Profiler::ObjectEnum(void){/*{{{*/ + + return ProfilerEnum; + +} +/*}}}*/ /*Profiler routines:*/ +IssmDouble Profiler::DeltaFlops(int inittag, int finaltag){/*{{{*/ + + /*Get initial flops*/ + _assert_(inittag>=0); + _assert_(inittag(this->flops[inittag])) _error_("Tag not set"); + IssmDouble init = this->flops[inittag]; + + /*Get final flops*/ + _assert_(finaltag>=0); + _assert_(finaltag(this->flops[finaltag])) _error_("Tag not set"); + IssmDouble final = this->flops[finaltag]; + + return final-init; +} +/*}}}*/ +IssmDouble Profiler::DeltaTime(int inittag, int finaltag){/*{{{*/ + + /*Get initial time*/ + _assert_(inittag>=0); + _assert_(inittag(this->time[inittag])) _error_("Tag not set"); + IssmDouble init = this->time[inittag]; + + /*Get final time*/ + _assert_(finaltag>=0); + _assert_(finaltag(this->time[finaltag])) _error_("Tag not set"); + IssmDouble final = this->time[finaltag]; + + #ifdef _HAVE_MPI_ + return final-init; + #else + return (final-init)/CLOCKS_PER_SEC; + #endif +} +/*}}}*/ +int Profiler::DeltaTimeModHour(int inittag, int finishtag){/*{{{*/ + + IssmDouble delta = this->DeltaTime(inittag,finishtag); + return int((reCast(delta))/3600); + +} +/*}}}*/ +int Profiler::DeltaTimeModMin(int inittag, int finishtag){/*{{{*/ + + IssmDouble delta = this->DeltaTime(inittag,finishtag); + return int(int(reCast(delta))%3600/60); +} +/*}}}*/ +int Profiler::DeltaTimeModSec(int inittag, int finishtag){/*{{{*/ + + IssmDouble delta = this->DeltaTime(inittag,finishtag); + return int(reCast(delta)%60); +} +/*}}}*/ +IssmDouble Profiler::Memory(int tag){/*{{{*/ + + /*Get initial flops*/ + _assert_(tag>=0); + _assert_(tag(this->flops[tag])) _error_("Tag not set"); + return this->memory[tag]; +} +/*}}}*/ void Profiler::Tag(int tagenum,bool dontmpisync){/*{{{*/ @@ -117,68 +183,2 @@ } /*}}}*/ -IssmDouble Profiler::DeltaTime(int inittag, int finaltag){/*{{{*/ - - /*Get initial time*/ - _assert_(inittag>=0); - _assert_(inittag(this->time[inittag])) _error_("Tag not set"); - IssmDouble init = this->time[inittag]; - - /*Get final time*/ - _assert_(finaltag>=0); - _assert_(finaltag(this->time[finaltag])) _error_("Tag not set"); - IssmDouble final = this->time[finaltag]; - - #ifdef _HAVE_MPI_ - return final-init; - #else - return (final-init)/CLOCKS_PER_SEC; - #endif -} -/*}}}*/ -IssmDouble Profiler::DeltaFlops(int inittag, int finaltag){/*{{{*/ - - /*Get initial flops*/ - _assert_(inittag>=0); - _assert_(inittag(this->flops[inittag])) _error_("Tag not set"); - IssmDouble init = this->flops[inittag]; - - /*Get final flops*/ - _assert_(finaltag>=0); - _assert_(finaltag(this->flops[finaltag])) _error_("Tag not set"); - IssmDouble final = this->flops[finaltag]; - - return final-init; -} -/*}}}*/ -int Profiler::DeltaTimeModHour(int inittag, int finishtag){/*{{{*/ - - IssmDouble delta = this->DeltaTime(inittag,finishtag); - return int((reCast(delta))/3600); - -} -/*}}}*/ -int Profiler::DeltaTimeModMin(int inittag, int finishtag){/*{{{*/ - - IssmDouble delta = this->DeltaTime(inittag,finishtag); - return int(int(reCast(delta))%3600/60); -} -/*}}}*/ -int Profiler::DeltaTimeModSec(int inittag, int finishtag){/*{{{*/ - - IssmDouble delta = this->DeltaTime(inittag,finishtag); - return int(reCast(delta)%60); -} -/*}}}*/ -IssmDouble Profiler::Memory(int tag){/*{{{*/ - - /*Get initial flops*/ - _assert_(tag>=0); - _assert_(tag(this->flops[tag])) _error_("Tag not set"); - return this->memory[tag]; -} -/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Profiler.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Profiler.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Profiler.h (revision 20810) @@ -27,7 +27,7 @@ public: - IssmDouble time[MAXIMUMSIZE]; IssmDouble flops[MAXIMUMSIZE]; IssmDouble memory[MAXIMUMSIZE]; + IssmDouble time[MAXIMUMSIZE]; /*Profiler constructors, destructors {{{*/ @@ -36,19 +36,19 @@ /*}}}*/ /*Object virtual functions definitions:{{{ */ + Object *copy(); + void DeepEcho(); void Echo(); - void DeepEcho(); int Id(); + void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction); int ObjectEnum(); - Object *copy(); - void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction); /*}}}*/ /*Profiler routines {{{*/ - void Tag(int tagenum,bool dontmpisync=false); - IssmDouble Memory(int tag); + IssmDouble DeltaFlops(int inittag, int finaltag); IssmDouble DeltaTime(int inittag, int finaltag); - IssmDouble DeltaFlops(int inittag, int finaltag); int DeltaTimeModHour(int inittag, int finaltag); int DeltaTimeModMin(int inittag, int finaltag); int DeltaTimeModSec(int inittag, int finaltag); + IssmDouble Memory(int tag); + void Tag(int tagenum,bool dontmpisync=false); /*}}}*/ }; Index: /issm/trunk-jpl/src/c/classes/RiftStruct.h =================================================================== --- /issm/trunk-jpl/src/c/classes/RiftStruct.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/RiftStruct.h (revision 20810) @@ -10,10 +10,10 @@ public: int numrifts; + int **riftspairs; + double **riftspenaltypairs; + int *riftsnumpairs; + int *riftsnumpenaltypairs; + int *riftsnumsegments; int **riftssegments; - int *riftsnumsegments; - int **riftspairs; - int *riftsnumpairs; - double **riftspenaltypairs; - int *riftsnumpenaltypairs; int *riftstips; double **state; Index: /issm/trunk-jpl/src/c/classes/Segment.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Segment.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Segment.h (revision 20810) @@ -46,4 +46,12 @@ /*Object virtual functions definitions:*/ + Object* copy() {/*{{{*/ + return new Segment(this->eid,this->x1,this->y1,this->x2,this->y2); + } + /*}}}*/ + void DeepEcho(void){/*{{{*/ + this->Echo(); + } + /*}}}*/ void Echo(void){/*{{{*/ @@ -55,9 +63,9 @@ } /*}}}*/ - void DeepEcho(void){/*{{{*/ - this->Echo(); - } + int Id(void){ return eid; }/*{{{*/ /*}}}*/ - int Id(void){ return eid; }/*{{{*/ + void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){/*{{{*/ + _error_("not implemented yet!"); + } /*}}}*/ int ObjectEnum(void){/*{{{*/ @@ -67,12 +75,4 @@ } /*}}}*/ - Object* copy() {/*{{{*/ - return new Segment(this->eid,this->x1,this->y1,this->x2,this->y2); - } - /*}}}*/ - void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){/*{{{*/ - _error_("not implemented yet!"); - } - /*}}}*/ }; Index: /issm/trunk-jpl/src/c/classes/Update.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Update.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Update.h (revision 20810) @@ -14,13 +14,13 @@ public: - virtual void InputUpdateFromVector(IssmDouble* vector, int name, int type)=0; + virtual void InputUpdateFromConstant(IssmDouble constant, int name)=0; + virtual void InputUpdateFromConstant(int constant, int name)=0; + virtual void InputUpdateFromConstant(bool constant, int name)=0; #ifdef _HAVE_DAKOTA_ virtual void InputUpdateFromMatrixDakota(IssmDouble* matrix, int rows, int ncols, int name, int type)=0; virtual void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type)=0; #endif - virtual void InputUpdateFromConstant(IssmDouble constant, int name)=0; - virtual void InputUpdateFromConstant(int constant, int name)=0; - virtual void InputUpdateFromConstant(bool constant, int name)=0; virtual void InputUpdateFromIoModel(int index, IoModel* iomodel)=0; + virtual void InputUpdateFromVector(IssmDouble* vector, int name, int type)=0; }; Index: /issm/trunk-jpl/src/c/classes/Vertex.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Vertex.cpp (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Vertex.cpp (revision 20810) @@ -62,4 +62,14 @@ /*Object virtual functions definitions:*/ +Object* Vertex::copy() {/*{{{*/ + + return new Vertex(*this); + +} +/*}}}*/ +void Vertex::DeepEcho(void){/*{{{*/ + this->Echo(); +} +/*}}}*/ void Vertex::Echo(void){/*{{{*/ @@ -78,21 +88,5 @@ } /*}}}*/ -void Vertex::DeepEcho(void){/*{{{*/ - this->Echo(); -} -/*}}}*/ int Vertex::Id(void){ return id; }/*{{{*/ -/*}}}*/ -int Vertex::ObjectEnum(void){/*{{{*/ - - return VertexEnum; - -} -/*}}}*/ -Object* Vertex::copy() {/*{{{*/ - - return new Vertex(*this); - -} /*}}}*/ void Vertex::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/ @@ -112,8 +106,42 @@ } /*}}}*/ +int Vertex::ObjectEnum(void){/*{{{*/ + + return VertexEnum; + +} +/*}}}*/ /*Vertex management: */ int Vertex::Connectivity(void){return connectivity;}/*{{{*/ /*}}}*/ +void Vertex::DistributePids(int* ppidcount){/*{{{*/ + + /*retrieve current pid*/ + int pidcount=*ppidcount; + + /*This vertex is a clone! Don't distribute pids, it will get them from another cpu!*/ + if(this->clone) return; + + /*This vertex should distribute its pid*/ + this->pid=pidcount; + pidcount++; + + /*Assign output pointers: */ + *ppidcount=pidcount; +} +/*}}}*/ +IssmDouble Vertex::GetLatitude(){/*{{{*/ + return this->latitute; +} +/*}}}*/ +IssmDouble Vertex::GetLongitude(){/*{{{*/ + return this->longitude; +} +/*}}}*/ +IssmDouble Vertex::GetRadius(){/*{{{*/ + return this->R; +} +/*}}}*/ IssmDouble Vertex::GetX(){/*{{{*/ return this->x; @@ -128,19 +156,70 @@ } /*}}}*/ -IssmDouble Vertex::GetLatitude(){/*{{{*/ - return this->latitute; -} -/*}}}*/ -IssmDouble Vertex::GetLongitude(){/*{{{*/ - return this->longitude; -} -/*}}}*/ -IssmDouble Vertex::GetRadius(){/*{{{*/ - return this->R; +void Vertex::OffsetPids(int pidcount){/*{{{*/ + + /*This vertex is a clone, don't offset the pids*/ + if(this->clone) return; + + /*This vertex should offset his pid, go ahead: */ + this->pid+=pidcount; +} +/*}}}*/ +int Vertex::Pid(void){ return pid; }/*{{{*/ +/*}}}*/ +void Vertex::SetClone(int* minranks){/*{{{*/ + + int my_rank; + + /*recover my_rank:*/ + my_rank=IssmComm::GetRank(); + + if (minranks[this->sid]==my_rank){ + this->clone=false; + } + else{ + /*!there is a cpu with lower rank that has the same vertex, + therefore, I am a clone*/ + this->clone=true; + } + +} +/*}}}*/ +void Vertex::ShowTruePids(int* truepids){/*{{{*/ + + /*Are we a clone? : */ + if(this->clone)return; + + /*Ok, we are not a clone, just plug our pid into truepids: */ + truepids[this->sid]=this->pid; } /*}}}*/ int Vertex::Sid(void){ return sid; }/*{{{*/ /*}}}*/ -int Vertex::Pid(void){ return pid; }/*{{{*/ +void Vertex::ToXYZ(Matrix* matrix){/*{{{*/ + + IssmDouble xyz[3]; + int indices[3]; + + if (this->clone==true) return; + + xyz[0]=x; + xyz[1]=y; + xyz[2]=z; + indices[0]=0; + indices[1]=1; + indices[2]=2; + + matrix->SetValues(1,&sid,3,&indices[0],&xyz[0],INS_VAL); +} +/*}}}*/ +void Vertex::UpdateClonePids(int* alltruepids){/*{{{*/ + + /*If we are not a clone, don't update, we already have pids: */ + if(!this->clone)return; + + /*Ok, we are a clone node, but we did not create the pid for this vertex + * Therefore, our pid is garbage right now. Go pick it up in the alltruepids: */ + this->pid=alltruepids[this->sid]; +} /*}}}*/ void Vertex::UpdatePosition(Vector* vx,Vector* vy,Vector* vz,Parameters* parameters,IssmDouble* surface,IssmDouble* bed){/*{{{*/ @@ -179,83 +258,4 @@ } /*}}}*/ -void Vertex::DistributePids(int* ppidcount){/*{{{*/ - - /*retrieve current pid*/ - int pidcount=*ppidcount; - - /*This vertex is a clone! Don't distribute pids, it will get them from another cpu!*/ - if(this->clone) return; - - /*This vertex should distribute its pid*/ - this->pid=pidcount; - pidcount++; - - /*Assign output pointers: */ - *ppidcount=pidcount; -} -/*}}}*/ -void Vertex::OffsetPids(int pidcount){/*{{{*/ - - /*This vertex is a clone, don't offset the pids*/ - if(this->clone) return; - - /*This vertex should offset his pid, go ahead: */ - this->pid+=pidcount; -} -/*}}}*/ -void Vertex::ShowTruePids(int* truepids){/*{{{*/ - - /*Are we a clone? : */ - if(this->clone)return; - - /*Ok, we are not a clone, just plug our pid into truepids: */ - truepids[this->sid]=this->pid; -} -/*}}}*/ -void Vertex::UpdateClonePids(int* alltruepids){/*{{{*/ - - /*If we are not a clone, don't update, we already have pids: */ - if(!this->clone)return; - - /*Ok, we are a clone node, but we did not create the pid for this vertex - * Therefore, our pid is garbage right now. Go pick it up in the alltruepids: */ - this->pid=alltruepids[this->sid]; -} -/*}}}*/ -void Vertex::SetClone(int* minranks){/*{{{*/ - - int my_rank; - - /*recover my_rank:*/ - my_rank=IssmComm::GetRank(); - - if (minranks[this->sid]==my_rank){ - this->clone=false; - } - else{ - /*!there is a cpu with lower rank that has the same vertex, - therefore, I am a clone*/ - this->clone=true; - } - -} -/*}}}*/ -void Vertex::ToXYZ(Matrix* matrix){/*{{{*/ - - IssmDouble xyz[3]; - int indices[3]; - - if (this->clone==true) return; - - xyz[0]=x; - xyz[1]=y; - xyz[2]=z; - indices[0]=0; - indices[1]=1; - indices[2]=2; - - matrix->SetValues(1,&sid,3,&indices[0],&xyz[0],INS_VAL); -} -/*}}}*/ void Vertex::VertexCoordinates(Vector* vx,Vector* vy,Vector* vz, bool spherical){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Vertex.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Vertex.h (revision 20809) +++ /issm/trunk-jpl/src/c/classes/Vertex.h (revision 20810) @@ -49,20 +49,20 @@ /*}}}*/ /*Vertex management:*/ - int Pid(void); - int Sid(void); int Connectivity(void); + void DistributePids(int* ppidcount); + IssmDouble GetLatitude(void); + IssmDouble GetLongitude(void); + IssmDouble GetRadius(void); IssmDouble GetX(void); IssmDouble GetY(void); IssmDouble GetZ(void); - IssmDouble GetLatitude(void); - IssmDouble GetLongitude(void); - IssmDouble GetRadius(void); + void OffsetPids(int pidcount); + int Pid(void); + void SetClone(int* minranks); + void ShowTruePids(int* borderpids); + int Sid(void); + void ToXYZ(Matrix* matrix); + void UpdateClonePids(int* allborderpids); void UpdatePosition(Vector* vx,Vector* vy,Vector* vz,Parameters* parameters,IssmDouble* thickness,IssmDouble* bed); - void DistributePids(int* ppidcount); - void OffsetPids(int pidcount); - void ShowTruePids(int* borderpids); - void UpdateClonePids(int* allborderpids); - void SetClone(int* minranks); - void ToXYZ(Matrix* matrix); void VertexCoordinates(Vector* vx,Vector* vy,Vector* vz,bool spherical=false); };