Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 15297) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 15298) @@ -87,6 +87,6 @@ ./classes/Inputs/Input.h\ ./classes/Inputs/InputLocal.h\ - ./classes/Inputs/TriaP1Input.h\ - ./classes/Inputs/TriaP1Input.cpp\ + ./classes/Inputs/TriaInput.h\ + ./classes/Inputs/TriaInput.cpp\ ./classes/Inputs/BoolInput.h\ ./classes/Inputs/BoolInput.cpp\ Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 15297) +++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 15298) @@ -1257,2 +1257,13 @@ } /*}}}*/ +/*FUNCTION PentaRef::NumberofNodes{{{*/ +int PentaRef::NumberofNodes(void){ + + switch(this->element_type){ + case P1Enum: return NUMNODESP1; + default: _error_("Element type "<element_type)<<" not supported yet"); + } + + return -1; +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h (revision 15297) +++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h (revision 15298) @@ -59,4 +59,5 @@ void GetInputDerivativeValue(IssmDouble* pvalues, IssmDouble* plist,IssmDouble* xyz_list, GaussTria* gauss){_error_("only PentaGauss are supported");}; + int NumberofNodes(void); }; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15297) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15298) @@ -887,10 +887,10 @@ /*Add Stress tensor components into inputs*/ - this->inputs->AddInput(new TriaP1Input(StressTensorxxEnum,&sigma_xx[0])); - this->inputs->AddInput(new TriaP1Input(StressTensorxyEnum,&sigma_xy[0])); - this->inputs->AddInput(new TriaP1Input(StressTensorxzEnum,&sigma_xz[0])); - this->inputs->AddInput(new TriaP1Input(StressTensoryyEnum,&sigma_yy[0])); - this->inputs->AddInput(new TriaP1Input(StressTensoryzEnum,&sigma_yz[0])); - this->inputs->AddInput(new TriaP1Input(StressTensorzzEnum,&sigma_zz[0])); + this->inputs->AddInput(new TriaInput(StressTensorxxEnum,&sigma_xx[0])); + this->inputs->AddInput(new TriaInput(StressTensorxyEnum,&sigma_xy[0])); + this->inputs->AddInput(new TriaInput(StressTensorxzEnum,&sigma_xz[0])); + this->inputs->AddInput(new TriaInput(StressTensoryyEnum,&sigma_yy[0])); + this->inputs->AddInput(new TriaInput(StressTensoryzEnum,&sigma_yz[0])); + this->inputs->AddInput(new TriaInput(StressTensorzzEnum,&sigma_zz[0])); /*Clean up and return*/ @@ -1028,9 +1028,9 @@ for (int imonth=0;imonth<12;imonth++) { for(i=0;iAddTimeInput(newmonthinput1,time+imonth/12.*yts); for(i=0;iAddTimeInput(newmonthinput2,time+imonth/12.*yts); } @@ -1394,4 +1394,45 @@ for (int iv=0;ivGetInput(enumtype); + int numnodes = this->NumberofNodes(); + + /* Start looping on the number of vertices: */ + if (input){ + GaussTria* gauss=new GaussTria(); + for (int iv=0;ivNumberofNodes();iv++){ + gauss->GaussNode(numnodes,iv); + input->GetInputValue(&pvalue[iv],gauss); + } + delete gauss; + } + else{ + for (int iv=0;ivGetInput(enumtype); + if (!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element"); + int numnodes = this->NumberofNodes(); + + /* Start looping on the number of vertices: */ + GaussTria* gauss=new GaussTria(); + for (int iv=0;ivNumberofNodes();iv++){ + gauss->GaussNode(numnodes,iv); + input->GetInputValue(&pvalue[iv],gauss); + } + delete gauss; } /*}}}*/ @@ -1687,5 +1728,5 @@ for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts; for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts; - this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); } break; @@ -1695,5 +1736,5 @@ for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts; for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts; - this->inputs->AddInput(new ControlInput(VxEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(VxEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); } break; @@ -1703,5 +1744,5 @@ for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts; for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts; - this->inputs->AddInput(new ControlInput(VyEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(VyEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); } break; @@ -1711,5 +1752,5 @@ for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]; for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]; - this->inputs->AddInput(new ControlInput(ThicknessEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(ThicknessEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); } break; @@ -1719,5 +1760,5 @@ for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]; for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]; - this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); } break; @@ -1739,5 +1780,5 @@ for(i=0;iData(InversionCostFunctionsCoefficientsEnum)[(tria_vertex_ids[j]-1)*num_cm_responses+i]; - datasetinput->inputs->AddObject(new TriaP1Input(InversionCostFunctionsCoefficientsEnum,nodeinputs)); + datasetinput->inputs->AddObject(new TriaInput(InversionCostFunctionsCoefficientsEnum,nodeinputs)); } @@ -1825,5 +1866,5 @@ /*Add input to the element: */ - this->inputs->AddInput(new TriaP1Input(enum_type,values)); + this->inputs->AddInput(new TriaInput(enum_type,values)); /*Free ressources:*/ @@ -1890,7 +1931,7 @@ /*Add input to the element: */ - this->inputs->AddInput(new TriaP1Input(ThicknessEnum,newthickness)); - this->inputs->AddInput(new TriaP1Input(SurfaceEnum,newsurface)); - this->inputs->AddInput(new TriaP1Input(BedEnum,newbed)); + this->inputs->AddInput(new TriaInput(ThicknessEnum,newthickness)); + this->inputs->AddInput(new TriaInput(SurfaceEnum,newsurface)); + this->inputs->AddInput(new TriaInput(BedEnum,newbed)); /*Free ressources:*/ @@ -1915,8 +1956,8 @@ /*update input*/ if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum || name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){ - material->inputs->AddInput(new TriaP1Input(name,values)); + material->inputs->AddInput(new TriaInput(name,values)); } else{ - this->inputs->AddInput(new TriaP1Input(name,values)); + this->inputs->AddInput(new TriaInput(name,values)); } return; @@ -1928,8 +1969,8 @@ /*update input*/ if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum || name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){ - material->inputs->AddInput(new TriaP1Input(name,values)); + material->inputs->AddInput(new TriaInput(name,values)); } else{ - this->inputs->AddInput(new TriaP1Input(name,values)); + this->inputs->AddInput(new TriaInput(name,values)); } return; @@ -1945,5 +1986,5 @@ } /*Add input to the element: */ - this->inputs->AddInput(new TriaP1Input(name,values)); + this->inputs->AddInput(new TriaInput(name,values)); /*Free ressources:*/ @@ -1957,5 +1998,5 @@ } /*Add input to the element: */ - this->inputs->AddInput(new TriaP1Input(name,values)); + this->inputs->AddInput(new TriaInput(name,values)); /*Free ressources:*/ @@ -2033,5 +2074,5 @@ /*create static input: */ - this->inputs->AddInput(new TriaP1Input(vector_enum,nodeinputs)); + this->inputs->AddInput(new TriaInput(vector_enum,nodeinputs)); } else if(M==numberofvertices+1){ @@ -2049,5 +2090,5 @@ if(t==0) transientinput=new TransientInput(vector_enum); - transientinput->AddTimeInput(new TriaP1Input(vector_enum,nodeinputs),time); + transientinput->AddTimeInput(new TriaInput(vector_enum,nodeinputs),time); } this->inputs->AddInput(transientinput); @@ -2316,16 +2357,16 @@ if(!this->IsFloating() && floatingelement==true){ for(i=0;iinputs->AddInput(new TriaP1Input(BasalforcingsMeltingRateEnum,&melting[0])); + this->inputs->AddInput(new TriaInput(BasalforcingsMeltingRateEnum,&melting[0])); } /*Update inputs*/ this->inputs->AddInput(new BoolInput(MaskElementonfloatingiceEnum,floatingelement)); - this->inputs->AddInput(new TriaP1Input(SurfaceEnum,&s[0])); - this->inputs->AddInput(new TriaP1Input(BedEnum,&b[0])); + this->inputs->AddInput(new TriaInput(SurfaceEnum,&s[0])); + this->inputs->AddInput(new TriaInput(BedEnum,&b[0])); /*Recalculate phi*/ if(migration_style==SubelementMigrationEnum || migration_style==SubelementMigration2Enum){ for(i=0;iinputs->AddInput(new TriaP1Input(GLlevelsetEnum,&phi[0])); + this->inputs->AddInput(new TriaInput(GLlevelsetEnum,&phi[0])); } } @@ -2486,5 +2527,5 @@ /*Update inputs*/ - this->inputs->AddInput(new TriaP1Input(SurfaceforcingsMassBalanceEnum,&agd[0])); + this->inputs->AddInput(new TriaInput(SurfaceforcingsMassBalanceEnum,&agd[0])); // this->inputs->AddInput(new TriaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0])); @@ -2576,5 +2617,5 @@ } //end of the loop over the vertices /*Update inputs*/ - this->inputs->AddInput(new TriaP1Input(SurfaceforcingsMassBalanceEnum,&smb[0])); + this->inputs->AddInput(new TriaInput(SurfaceforcingsMassBalanceEnum,&smb[0])); } /*}}}*/ @@ -2760,22 +2801,22 @@ if(!iomodel->Data(VxEnum)){ for(i=0;i<3;i++)nodeinputs[i]=0; - this->inputs->AddInput(new TriaP1Input(VxEnum,nodeinputs)); - if(dakota_analysis) this->inputs->AddInput(new TriaP1Input(QmuVxEnum,nodeinputs)); + this->inputs->AddInput(new TriaInput(VxEnum,nodeinputs)); + if(dakota_analysis) this->inputs->AddInput(new TriaInput(QmuVxEnum,nodeinputs)); } if(!iomodel->Data(VyEnum)){ for(i=0;i<3;i++)nodeinputs[i]=0; - this->inputs->AddInput(new TriaP1Input(VyEnum,nodeinputs)); - if(dakota_analysis) this->inputs->AddInput(new TriaP1Input(QmuVyEnum,nodeinputs)); + this->inputs->AddInput(new TriaInput(VyEnum,nodeinputs)); + if(dakota_analysis) this->inputs->AddInput(new TriaInput(QmuVyEnum,nodeinputs)); } if(!iomodel->Data(VzEnum)){ for(i=0;i<3;i++)nodeinputs[i]=0; - this->inputs->AddInput(new TriaP1Input(VzEnum,nodeinputs)); - if(dakota_analysis) this->inputs->AddInput(new TriaP1Input(QmuVzEnum,nodeinputs)); + this->inputs->AddInput(new TriaInput(VzEnum,nodeinputs)); + if(dakota_analysis) this->inputs->AddInput(new TriaInput(QmuVzEnum,nodeinputs)); } if(!iomodel->Data(PressureEnum)){ for(i=0;i<3;i++)nodeinputs[i]=0; if(dakota_analysis){ - this->inputs->AddInput(new TriaP1Input(PressureEnum,nodeinputs)); - this->inputs->AddInput(new TriaP1Input(QmuPressureEnum,nodeinputs)); + this->inputs->AddInput(new TriaInput(PressureEnum,nodeinputs)); + this->inputs->AddInput(new TriaInput(QmuPressureEnum,nodeinputs)); } } @@ -3236,7 +3277,4 @@ ElementMatrix* Tria::CreateKMatrixDiagnosticMacAyealViscous(void){ - /*Constants*/ - const int numdof=NDOF2*NUMVERTICES; - /*Intermediaries*/ int i,j; @@ -3245,12 +3283,15 @@ IssmDouble viscosity_overshoot,thickness,Jdet; IssmDouble epsilon[3],oldepsilon[3]; /* epsilon=[exx,eyy,exy]; */ - IssmDouble B[3][numdof]; - IssmDouble Bprime[3][numdof]; - IssmDouble D[3][3] = {0.0}; IssmDouble D_scalar; - GaussTria *gauss = NULL; - - /*Initialize Element matrix*/ - ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum); + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); + int numdof = numnodes*NDOF2; + + /*Initialize Element matrix, vectors and Gaussian points*/ + ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,MacAyealApproximationEnum); + IssmDouble* B = xNew(3*numdof); + IssmDouble* Bprime = xNew(3*numdof); + IssmDouble* D = xNewZeroInit(3*3); /*Retrieve all inputs and parameters*/ @@ -3264,5 +3305,5 @@ /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(2); + GaussTria* gauss = new GaussTria(2); for(int ig=gauss->begin();igend();ig++){ @@ -3270,6 +3311,6 @@ GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss); - GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss); + GetBMacAyeal(&B[0], &xyz_list[0][0], gauss); + GetBprimeMacAyeal(&Bprime[0], &xyz_list[0][0], gauss); this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); @@ -3280,10 +3321,10 @@ newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); - D_scalar=2*newviscosity*thickness*gauss->weight*Jdet; - for (i=0;i<3;i++) D[i][i]=D_scalar; - - TripleMultiply(&B[0][0],3,numdof,1, - &D[0][0],3,3,0, - &Bprime[0][0],3,numdof,0, + D_scalar=2.*newviscosity*thickness*gauss->weight*Jdet; + for(i=0;i<3;i++) D[i*3+i]=D_scalar; + + TripleMultiply(B,3,numdof,1, + D,3,3,0, + Bprime,3,numdof,0, &Ke->values[0],1); } @@ -3294,4 +3335,7 @@ /*Clean up and return*/ delete gauss; + xDelete(B); + xDelete(D); + xDelete(Bprime); return Ke; } @@ -3299,7 +3343,4 @@ /*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealFriction {{{*/ ElementMatrix* Tria::CreateKMatrixDiagnosticMacAyealFriction(void){ - - /*Constants*/ - const int numdof=NDOF2*NUMVERTICES; /*Intermediaries*/ @@ -3312,8 +3353,5 @@ IssmDouble fraction1,fraction2; IssmDouble phi=1.0; - IssmDouble L[2][numdof]; - IssmDouble DL[2][2] = {{ 0,0 },{0,0}}; - IssmDouble DL_scalar; - IssmDouble slope[2] = {0.0,0.0}; + IssmDouble D_scalar; IssmDouble gllevelset; IssmDouble xyz_list[NUMVERTICES][3]; @@ -3321,7 +3359,15 @@ GaussTria *gauss = NULL; - /*Initialize Element matrix and return if necessary*/ + /*Return if element is inactive*/ if(IsFloating()) return NULL; - ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum); + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); + int numdof = numnodes*NDOF2; + + /*Initialize Element matrix and vectors*/ + ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,MacAyealApproximationEnum); + IssmDouble* B = xNew(2*numdof); + IssmDouble* D = xNewZeroInit(2*2); /*Retrieve all inputs and parameters*/ @@ -3361,12 +3407,12 @@ } - GetL(&L[0][0], &xyz_list[0][0], gauss,NDOF2); + GetL(&B[0], &xyz_list[0][0], gauss,NDOF2); GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - DL_scalar=alpha2*gauss->weight*Jdet; - for (i=0;i<2;i++) DL[i][i]=DL_scalar; - - TripleMultiply( &L[0][0],2,numdof,1, - &DL[0][0],2,2,0, - &L[0][0],2,numdof,0, + D_scalar=alpha2*gauss->weight*Jdet; + for(i=0;i<2;i++) D[i*2+i]=D_scalar; + + TripleMultiply(B,2,numdof,1, + D,2,2,0, + B,2,numdof,0, &Ke->values[0],1); } @@ -3378,4 +3424,6 @@ delete gauss; delete friction; + xDelete(B); + xDelete(D); return Ke; } @@ -3404,7 +3452,4 @@ /*FUNCTION Tria::CreatePVectorDiagnosticMacAyeal {{{*/ ElementVector* Tria::CreatePVectorDiagnosticMacAyeal(){ - - /*Constants*/ - const int numdof=NDOF2*NUMVERTICES; /*Intermediaries */ @@ -3414,10 +3459,13 @@ IssmDouble xyz_list[NUMVERTICES][3]; IssmDouble slope[2]; - IssmDouble basis[3]; - IssmDouble pe_g_gaussian[numdof]; - GaussTria* gauss=NULL; - - /*Initialize Element vector*/ + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); + int numdof = numnodes*NDOF2; + + /*Initialize Element vector and vectors*/ ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum); + GaussTria* gauss = new GaussTria(2); + IssmDouble* basis = xNew(numnodes); /*Retrieve all inputs and parameters*/ @@ -3428,5 +3476,4 @@ /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(2); for(int ig=gauss->begin();igend();ig++){ @@ -3440,6 +3487,6 @@ driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG()*thickness; - /*Build pe_g_gaussian vector: */ - for (i=0;ivalues[i*NDOF2+j]+=-driving_stress_baseline*slope[j]*Jdet*gauss->weight*basis[i]; @@ -3453,4 +3500,5 @@ /*Clean up and return*/ delete gauss; + xDelete(basis); return pe; } @@ -3508,7 +3556,4 @@ /*FUNCTION Tria::CreateJacobianDiagnosticMacayeal{{{*/ ElementMatrix* Tria::CreateJacobianDiagnosticMacayeal(void){ - - /*Constants*/ - const int numdof=NDOF2*NUMVERTICES; /*Intermediaries */ @@ -3521,10 +3566,13 @@ IssmDouble epsilon[3];/* epsilon=[exx,eyy,exy];*/ IssmDouble eps1[2],eps2[2]; - IssmDouble phi[NUMVERTICES]; - IssmDouble dphi[2][NUMVERTICES]; - GaussTria *gauss=NULL; - - /*Initialize Jacobian with regular MacAyeal (first part of the Gateau derivative)*/ - ElementMatrix* Ke=CreateKMatrixDiagnosticMacAyeal(); + GaussTria* gauss = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); + int numdof = numnodes*NDOF2; + + /*Initialize Element matrix, vectors and Gaussian points*/ + ElementMatrix* Ke=CreateKMatrixDiagnosticMacAyeal(); //Initialize Jacobian with regular MacAyeal (first part of the Gateau derivative) + IssmDouble* dphi = xNew(2*numnodes); /*Retrieve all inputs and parameters*/ @@ -3541,5 +3589,5 @@ GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctionsDerivatives(&dphi[0][0],&xyz_list[0][0],gauss); + GetNodalFunctionsDerivatives(dphi,&xyz_list[0][0],gauss); thickness_input->GetInputValue(&thickness, gauss); @@ -3549,15 +3597,15 @@ eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1]; - for(i=0;i<3;i++){ - for(j=0;j<3;j++){ - eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]; - eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]; - eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]; - eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]; - - Ke->values[6*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii; - Ke->values[6*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii; - Ke->values[6*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii; - Ke->values[6*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii; + for(i=0;ivalues[2*numnodes*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii; + Ke->values[2*numnodes*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii; + Ke->values[2*numnodes*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii; + Ke->values[2*numnodes*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii; } } @@ -3568,4 +3616,5 @@ /*Clean up and return*/ + xDelete(dphi); delete gauss; return Ke; @@ -3575,14 +3624,15 @@ void Tria::GetSolutionFromInputsDiagnosticHoriz(Vector* solution){ - const int numdof=NDOF2*NUMVERTICES; - - int i; - int* doflist=NULL; - IssmDouble vx,vy; - IssmDouble values[numdof]; - GaussTria* gauss=NULL; - - /*Get dof list: */ + IssmDouble vx,vy; + int* doflist = NULL; + GaussTria* gauss = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); + int numdof = numnodes*NDOF2; + + /*Fetch dof list and allocate solution vector*/ GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + IssmDouble* values = xNew(numdof); /*Get inputs*/ @@ -3593,7 +3643,7 @@ /*P1 element only for now*/ gauss=new GaussTria(); - for(i=0;iGaussVertex(i); + for(int i=0;iGaussNode(numnodes,i); /*Recover vx and vy*/ @@ -3609,4 +3659,5 @@ delete gauss; xDelete(doflist); + xDelete(values); } /*}}}*/ @@ -3653,19 +3704,21 @@ void Tria::InputUpdateFromSolutionDiagnosticHoriz(IssmDouble* solution){ - const int numdof=NDOF2*NUMVERTICES; - - int i; - int* doflist=NULL; - IssmDouble rho_ice,g; - IssmDouble values[numdof]; - IssmDouble vx[NUMVERTICES]; - IssmDouble vy[NUMVERTICES]; - IssmDouble vz[NUMVERTICES]; - IssmDouble vel[NUMVERTICES]; - IssmDouble pressure[NUMVERTICES]; - IssmDouble thickness[NUMVERTICES]; - - /*Get dof list: */ + int i; + IssmDouble rho_ice,g; + int* doflist=NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); + int numdof = numnodes*NDOF2; + + /*Fetch dof list and allocate solution vectors*/ GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + IssmDouble* values = xNew(numdof); + IssmDouble* vx = xNew(numdof); + IssmDouble* vy = xNew(numdof); + IssmDouble* vz = xNew(numdof); + IssmDouble* vel = xNew(numdof); + IssmDouble* pressure = xNew(numdof); + IssmDouble* thickness = xNew(numdof); /*Use the dof list to index into the solution vector: */ @@ -3673,8 +3726,8 @@ /*Transform solution in Cartesian Space*/ - TransformSolutionCoord(&values[0],nodes,NUMVERTICES,XYEnum); + TransformSolutionCoord(&values[0],nodes,numnodes,XYEnum); /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ - for(i=0;iGetRhoIce(); g=matpar->GetG(); - GetInputListOnVertices(&thickness[0],ThicknessEnum); - for(i=0;iinputs->AddInput(new TriaP1Input(VxEnum,vx)); - this->inputs->AddInput(new TriaP1Input(VyEnum,vy)); - this->inputs->AddInput(new TriaP1Input(VelEnum,vel)); - this->inputs->AddInput(new TriaP1Input(PressureEnum,pressure)); + this->inputs->AddInput(new TriaInput(VxEnum,vx)); + this->inputs->AddInput(new TriaInput(VyEnum,vy)); + this->inputs->AddInput(new TriaInput(VelEnum,vel)); + this->inputs->AddInput(new TriaInput(PressureEnum,pressure)); /*Free ressources:*/ xDelete(doflist); - + xDelete(values); + xDelete(vx); + xDelete(vy); + xDelete(vz); + xDelete(vel); + xDelete(pressure); + xDelete(thickness); } /*}}}*/ @@ -3716,25 +3775,30 @@ void Tria::InputUpdateFromSolutionDiagnosticHutter(IssmDouble* solution){ - const int numdof=NDOF2*NUMVERTICES; - - int i; - int* doflist=NULL; - IssmDouble rho_ice,g; - IssmDouble values[numdof]; - IssmDouble vx[NUMVERTICES]; - IssmDouble vy[NUMVERTICES]; - IssmDouble vz[NUMVERTICES]; - IssmDouble vel[NUMVERTICES]; - IssmDouble pressure[NUMVERTICES]; - IssmDouble thickness[NUMVERTICES]; - - /*Get dof list: */ + int i; + IssmDouble rho_ice,g; + int* doflist=NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); + int numdof = numnodes*NDOF2; + + /*Fetch dof list and allocate solution vectors*/ GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + IssmDouble* values = xNew(numdof); + IssmDouble* vx = xNew(numdof); + IssmDouble* vy = xNew(numdof); + IssmDouble* vz = xNew(numdof); + IssmDouble* vel = xNew(numdof); + IssmDouble* pressure = xNew(numdof); + IssmDouble* thickness = xNew(numdof); /*Use the dof list to index into the solution vector: */ for(i=0;iGetRhoIce(); g=matpar->GetG(); - GetInputListOnVertices(&thickness[0],ThicknessEnum); - for(i=0;iinputs->AddInput(new TriaP1Input(VxEnum,vx)); - this->inputs->AddInput(new TriaP1Input(VyEnum,vy)); - this->inputs->AddInput(new TriaP1Input(VelEnum,vel)); - this->inputs->AddInput(new TriaP1Input(PressureEnum,pressure)); + this->inputs->AddInput(new TriaInput(VxEnum,vx)); + this->inputs->AddInput(new TriaInput(VyEnum,vy)); + this->inputs->AddInput(new TriaInput(VelEnum,vel)); + this->inputs->AddInput(new TriaInput(PressureEnum,pressure)); /*Free ressources:*/ xDelete(doflist); + xDelete(values); + xDelete(vx); + xDelete(vy); + xDelete(vz); + xDelete(vel); + xDelete(pressure); + xDelete(thickness); } /*}}}*/ @@ -3922,5 +3993,5 @@ GradientIndexing(&vertexpidlist[0],control_index); for(int i=0;iSetGradient(grad_input); @@ -5689,6 +5760,6 @@ /*Add vx and vy as inputs to the tria element: */ - this->inputs->AddInput(new TriaP1Input(AdjointxEnum,lambdax)); - this->inputs->AddInput(new TriaP1Input(AdjointyEnum,lambday)); + this->inputs->AddInput(new TriaInput(AdjointxEnum,lambdax)); + this->inputs->AddInput(new TriaInput(AdjointyEnum,lambday)); /*Free ressources:*/ @@ -5719,5 +5790,5 @@ /*Add vx and vy as inputs to the tria element: */ - this->inputs->AddInput(new TriaP1Input(AdjointEnum,lambda)); + this->inputs->AddInput(new TriaInput(AdjointEnum,lambda)); /*Free ressources:*/ @@ -5771,5 +5842,5 @@ values[i]=vector[vertexpidlist[i]]; } - new_input = new TriaP1Input(control_enum,values); + new_input = new TriaInput(control_enum,values); if(control_enum==MaterialsRheologyBbarEnum || control_enum==MaterialsRheologyZbarEnum){ @@ -5865,6 +5936,6 @@ /*Add to inputs*/ - this->inputs->AddInput(new TriaP1Input(HydrologyWaterVxEnum,vx)); - this->inputs->AddInput(new TriaP1Input(HydrologyWaterVyEnum,vy)); + this->inputs->AddInput(new TriaInput(HydrologyWaterVxEnum,vx)); + this->inputs->AddInput(new TriaInput(HydrologyWaterVyEnum,vy)); } /*}}}*/ @@ -6399,5 +6470,5 @@ /*Add input to the element: */ - this->inputs->AddInput(new TriaP1Input(WatercolumnEnum,values)); + this->inputs->AddInput(new TriaInput(WatercolumnEnum,values)); /*Free ressources:*/ @@ -6453,6 +6524,6 @@ /*Add input to the element: */ - this->inputs->AddInput(new TriaP1Input(SedimentHeadEnum,values)); - this->inputs->AddInput(new TriaP1Input(SedimentHeadResidualEnum,residual)); + this->inputs->AddInput(new TriaInput(SedimentHeadEnum,values)); + this->inputs->AddInput(new TriaInput(SedimentHeadResidualEnum,residual)); /*Free ressources:*/ @@ -6657,5 +6728,5 @@ case VertexEnum: - /*New TriaP1Input*/ + /*New TriaInput*/ IssmDouble values[3]; @@ -6718,11 +6789,11 @@ /*Add new inputs: */ - this->inputs->AddInput(new TriaP1Input(ThicknessEnum,thickness)); - this->inputs->AddInput(new TriaP1Input(BedEnum,bed)); - this->inputs->AddInput(new TriaP1Input(SurfaceEnum,surface)); + this->inputs->AddInput(new TriaInput(ThicknessEnum,thickness)); + this->inputs->AddInput(new TriaInput(BedEnum,bed)); + this->inputs->AddInput(new TriaInput(SurfaceEnum,surface)); break; default: - this->inputs->AddInput(new TriaP1Input(name,values)); + this->inputs->AddInput(new TriaInput(name,values)); } break; @@ -6776,5 +6847,5 @@ if(t==0) transientinput=new TransientInput(name); - transientinput->AddTimeInput(new TriaP1Input(name,values),time); + transientinput->AddTimeInput(new TriaInput(name,values),time); transientinput->Configure(parameters); } Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 15297) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 15298) @@ -206,4 +206,6 @@ void GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue); void GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue,int index); //TO BE REMOVED + void GetInputListOnNodes(IssmDouble* pvalue,int enumtype); + void GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue); void GetInputValue(IssmDouble* pvalue,Node* node,int enumtype); void GetStrainRate2d(IssmDouble* epsilon,IssmDouble* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input); Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 15297) +++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 15298) @@ -26,5 +26,4 @@ /*}}}*/ /*FUNCTION TriaRef::TriaRef(int* types,int nummodels){{{*/ - TriaRef::TriaRef(const int nummodels){ @@ -222,10 +221,10 @@ /*Build B': */ for (int i=0;icoord1; - basis[1]=gauss->coord2; - basis[2]=gauss->coord3; - -} -/*}}}*/ -/*FUNCTION TriaRef::GetNodalFunctionsP2{{{*/ -void TriaRef::GetNodalFunctionsP2(IssmDouble** pbasis,GaussTria* gauss){ - /*This routine returns the values of the nodal functions at the gaussian point.*/ - - /*Allocate*/ - IssmDouble* basis =xNew(NUMNODESP2); - - basis[0]=gauss->coord1*(2.*gauss->coord1-1.); - basis[1]=gauss->coord2*(2.*gauss->coord2-1.); - basis[2]=gauss->coord3*(2.*gauss->coord3-1.); - basis[3]=4.*gauss->coord3*gauss->coord2; - basis[4]=4.*gauss->coord3*gauss->coord1; - basis[5]=4.*gauss->coord1*gauss->coord2; - - /*Assign output pointer*/ - *pbasis = basis; - + switch(this->element_type){ + case P1Enum: + case P1DGEnum: + basis[0]=gauss->coord1; + basis[1]=gauss->coord2; + basis[2]=gauss->coord3; + return; + case P2Enum: + basis[0]=gauss->coord1*(2.*gauss->coord1-1.); + basis[1]=gauss->coord2*(2.*gauss->coord2-1.); + basis[2]=gauss->coord3*(2.*gauss->coord3-1.); + basis[3]=4.*gauss->coord3*gauss->coord2; + basis[4]=4.*gauss->coord3*gauss->coord1; + basis[5]=4.*gauss->coord1*gauss->coord2; + return; + default: + _error_("Element type "<element_type)<<" not supported yet"); + } } /*}}}*/ @@ -424,12 +418,31 @@ /*This routine returns the values of the nodal functions at the gaussian point.*/ - IssmDouble BasisFunctions[3]; - - GetNodalFunctions(&BasisFunctions[0],gauss); - _assert_(index1>=0 && index1<3); _assert_(index2>=0 && index2<3); - basis[0]=BasisFunctions[index1]; - basis[1]=BasisFunctions[index2]; + + /*Fetch number of nodes for this finite element*/ + int numnodes = this->NumberofNodes(); + + /*Get nodal functions*/ + IssmDouble* BasisFunctions=xNew(numnodes); + GetNodalFunctions(BasisFunctions,gauss); + + switch(this->element_type){ + case P1Enum: + case P1DGEnum: + basis[0]=BasisFunctions[index1]; + basis[1]=BasisFunctions[index2]; + return; + case P2Enum: + _assert_(index2element_type)<<" not supported yet"); + } + + xDelete(BasisFunctions); } /*}}}*/ @@ -439,10 +452,12 @@ /*This routine returns the values of the nodal functions derivatives (with respect to the * actual coordinate system): */ - int i; - IssmDouble dbasis_ref[NDOF2][NUMNODESP1]; IssmDouble Jinv[NDOF2][NDOF2]; - /*Get derivative values with respect to parametric coordinate system: */ - GetNodalFunctionsDerivativesReference(&dbasis_ref[0][0], gauss); + /*Fetch number of nodes for this finite element*/ + int numnodes = this->NumberofNodes(); + + /*Get nodal functions derivatives in reference triangle*/ + IssmDouble* dbasis_ref=xNew(2*numnodes); + GetNodalFunctionsDerivativesReference(dbasis_ref,gauss); /*Get Jacobian invert: */ @@ -450,11 +465,10 @@ /*Build dbasis: - * * [dhi/dx]= Jinv*[dhi/dr] * [dhi/dy] [dhi/ds] */ - for (i=0;i(NUMNODESP2*2); - - /*Nodal function 1*/ - dbasis[NUMNODESP2*0+0]=-2.*gauss->coord1 + 0.5; - dbasis[NUMNODESP2*0+1]=-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.; - /*Nodal function 2*/ - dbasis[NUMNODESP2*0+0]=+2.*gauss->coord2 + 0.5; - dbasis[NUMNODESP2*0+1]=-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.; - /*Nodal function 3*/ - dbasis[NUMNODESP2*0+0]=0.; - dbasis[NUMNODESP2*0+1]=+4.*SQRT3/3.*gauss->coord3 - SQRT3/3.; - /*Nodal function 4*/ - dbasis[NUMNODESP2*0+0]=+2.*gauss->coord3; - dbasis[NUMNODESP2*0+1]=+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3; - /*Nodal function 5*/ - dbasis[NUMNODESP2*0+0]=-2.*gauss->coord3; - dbasis[NUMNODESP2*0+1]=+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3; - /*Nodal function 6*/ - dbasis[NUMNODESP2*0+0]=2.*(gauss->coord1-gauss->coord2); - dbasis[NUMNODESP2*0+1]=-2.*SQRT3/3.*(gauss->coord1+gauss->coord2); - - /*Assign output pointer*/ - *pdbasis = dbasis; + switch(this->element_type){ + case P1Enum: case P1DGEnum: + /*Nodal function 1*/ + dbasis[NUMNODESP1*0+0]=-0.5; + dbasis[NUMNODESP1*1+0]=-1.0/(2.0*SQRT3); + /*Nodal function 2*/ + dbasis[NUMNODESP1*0+1]=0.5; + dbasis[NUMNODESP1*1+1]=-1.0/(2.0*SQRT3); + /*Nodal function 3*/ + dbasis[NUMNODESP1*0+2]=0; + dbasis[NUMNODESP1*1+2]=1.0/SQRT3; + return; + case P2Enum: + /*Nodal function 1*/ + dbasis[NUMNODESP2*0+0]=-2.*gauss->coord1 + 0.5; + dbasis[NUMNODESP2*1+0]=-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.; + /*Nodal function 2*/ + dbasis[NUMNODESP2*0+1]=+2.*gauss->coord2 + 0.5; + dbasis[NUMNODESP2*1+1]=-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.; + /*Nodal function 3*/ + dbasis[NUMNODESP2*0+2]=0.; + dbasis[NUMNODESP2*1+2]=+4.*SQRT3/3.*gauss->coord3 - SQRT3/3.; + /*Nodal function 4*/ + dbasis[NUMNODESP2*0+3]=+2.*gauss->coord3; + dbasis[NUMNODESP2*1+3]=+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3; + /*Nodal function 5*/ + dbasis[NUMNODESP2*0+4]=-2.*gauss->coord3; + dbasis[NUMNODESP2*1+4]=+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3; + /*Nodal function 6*/ + dbasis[NUMNODESP2*0+5]=2.*(gauss->coord1-gauss->coord2); + dbasis[NUMNODESP2*1+5]=-2.*SQRT3/3.*(gauss->coord1+gauss->coord2); + return; + default: + _error_("Element type "<element_type)<<" not supported yet"); + } + } /*}}}*/ @@ -519,16 +526,25 @@ * dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx * - * p is a vector of size 2x1 already allocated. - */ - - /*Nodal Derivatives*/ - IssmDouble dbasis[2][3]; //nodal derivative functions in actual coordinate system. - - /*Get dh1dh2dh3 in actual coordinate system: */ - GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list, gauss); + * p is a vector already allocated. + */ + + /*Output*/ + IssmDouble dpx=0.; + IssmDouble dpy=0.; + + /*Fetch number of nodes for this finite element*/ + int numnodes = this->NumberofNodes(); + + /*Get nodal functions derivatives*/ + IssmDouble* dbasis=xNew(2*numnodes); + GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss); + + /*Calculate parameter for this Gauss point*/ + for(int i=0;iNumberofNodes(); /*Get nodal functions*/ + IssmDouble* basis=xNew(numnodes); GetNodalFunctions(basis, gauss); - /*Get parameter*/ - *p=basis[0]*plist[0]+basis[1]*plist[1]+basis[2]*plist[2]; + /*Calculate parameter for this Gauss point*/ + for(int i=0;ielement_type)<<" not supported yet"); + default: _error_("Element type "<element_type)<<" not supported yet"); } Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h (revision 15297) +++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h (revision 15298) @@ -36,5 +36,4 @@ void GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,GaussTria* gauss); void GetNodalFunctions(IssmDouble* basis,GaussTria* gauss); - void GetNodalFunctionsP2(IssmDouble** pbasis,GaussTria* gauss); void GetSegmentNodalFunctions(IssmDouble* basis,GaussTria* gauss, int index1,int index2); void GetSegmentBFlux(IssmDouble* B,GaussTria* gauss, int index1,int index2); @@ -42,10 +41,8 @@ void GetNodalFunctionsDerivatives(IssmDouble* basis,IssmDouble* xyz_list, GaussTria* gauss); void GetNodalFunctionsDerivativesReference(IssmDouble* dl1dl3,GaussTria* gauss); - void GetNodalFunctionsDerivativesP2Reference(IssmDouble** pdbasis,GaussTria* gauss); void GetInputValue(IssmDouble* pp, IssmDouble* plist, GaussTria* gauss); void GetInputDerivativeValue(IssmDouble* pp, IssmDouble* plist,IssmDouble* xyz_list, GaussTria* gauss); int NumberofNodes(void); - }; #endif Index: /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.cpp (revision 15297) +++ /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.cpp (revision 15298) @@ -30,9 +30,9 @@ switch(enum_input){ - case TriaP1InputEnum: - values =new TriaP1Input(enum_type,pvalues); - savedvalues=new TriaP1Input(enum_type,pvalues); - minvalues =new TriaP1Input(enum_type,pmin); - maxvalues =new TriaP1Input(enum_type,pmax); + case TriaInputEnum: + values =new TriaInput(enum_type,pvalues); + savedvalues=new TriaInput(enum_type,pvalues); + minvalues =new TriaInput(enum_type,pmin); + maxvalues =new TriaInput(enum_type,pmax); break; case PentaP1InputEnum: Index: /issm/trunk-jpl/src/c/classes/Inputs/PentaP1Input.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/PentaP1Input.cpp (revision 15297) +++ /issm/trunk-jpl/src/c/classes/Inputs/PentaP1Input.cpp (revision 15298) @@ -86,5 +86,5 @@ /*output*/ - TriaP1Input* outinput=NULL; + TriaInput* outinput=NULL; IssmDouble newvalues[3]; @@ -100,5 +100,5 @@ /*Create new Tria input*/ - outinput=new TriaP1Input(this->enum_type,&newvalues[0]); + outinput=new TriaInput(this->enum_type,&newvalues[0]); /*Assign output*/ Index: /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp (revision 15298) +++ /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp (revision 15298) @@ -0,0 +1,444 @@ +/*!\file TriaInput.c + * \brief: implementation of the TriaInput object + */ + +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../classes.h" +#include "../../shared/shared.h" + +/*TriaInput constructors and destructor*/ +/*FUNCTION TriaInput::TriaInput(){{{*/ +TriaInput::TriaInput(){ + values = NULL; +} +/*}}}*/ +/*FUNCTION TriaInput::TriaInput(int in_enum_type,IssmDouble* invalues,element_type_in){{{*/ +TriaInput::TriaInput(int in_enum_type,IssmDouble* in_values,int element_type_in) + :TriaRef(1) +{ + + /*Set TriaRef*/ + this->SetElementType(element_type_in,0); + this->element_type=element_type_in; + + /*Set Enum*/ + enum_type=in_enum_type; + + /*Set values*/ + this->values=xNew(this->NumberofNodes()); + for(int i=0;iNumberofNodes();i++) values[i]=in_values[i]; +} +/*}}}*/ +/*FUNCTION TriaInput::~TriaInput(){{{*/ +TriaInput::~TriaInput(){ + xDelete(this->values); + return; +} +/*}}}*/ + +/*Object virtual functions definitions:*/ +/*FUNCTION TriaInput::Echo {{{*/ +void TriaInput::Echo(void){ + this->DeepEcho(); +} +/*}}}*/ +/*FUNCTION TriaInput::DeepEcho{{{*/ +void TriaInput::DeepEcho(void){ + + _printf_("TriaInput:\n"); + _printf_(" enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")\n"); + _printf_(" values: ["); + for(int i=0;iNumberofNodes();i++) _printf_(" "<values[i]); + _printf_("]\n"); +} +/*}}}*/ +/*FUNCTION TriaInput::Id{{{*/ +int TriaInput::Id(void){ return -1; } +/*}}}*/ +/*FUNCTION TriaInput::ObjectEnum{{{*/ +int TriaInput::ObjectEnum(void){ + + return TriaInputEnum; + +} +/*}}}*/ +/*FUNCTION TriaInput::copy{{{*/ +Object* TriaInput::copy() { + + return new TriaInput(this->enum_type,this->values,this->element_type); + +} +/*}}}*/ + +/*TriaInput management*/ +/*FUNCTION TriaInput::InstanceEnum{{{*/ +int TriaInput::InstanceEnum(void){ + + return this->enum_type; + +} +/*}}}*/ +/*FUNCTION TriaInput::SpawnTriaInput{{{*/ +Input* TriaInput::SpawnTriaInput(int* indices){ + + /*output*/ + TriaInput* outinput=NULL; + + /*Create new Tria input (copy of current input)*/ + outinput=new TriaInput(this->enum_type,&this->values[0],this->element_type); + + /*Assign output*/ + return outinput; + +} +/*}}}*/ +/*FUNCTION TriaInput::SpawnResult{{{*/ +ElementResult* TriaInput::SpawnResult(int step, IssmDouble time){ + + return new TriaP1ElementResult(this->enum_type,this->values,step,time); + +} +/*}}}*/ + +/*Object functions*/ +/*FUNCTION TriaInput::GetInputValue(IssmDouble* pvalue,GaussTria* gauss){{{*/ +void TriaInput::GetInputValue(IssmDouble* pvalue,GaussTria* gauss){ + + /*Call TriaRef function*/ + TriaRef::GetInputValue(pvalue,&values[0],gauss); + +} +/*}}}*/ +/*FUNCTION TriaInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, GaussTria* gauss){{{*/ +void TriaInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, GaussTria* gauss){ + + /*Call TriaRef function*/ + TriaRef::GetInputDerivativeValue(p,&values[0],xyz_list,gauss); +} +/*}}}*/ +/*FUNCTION TriaInput::GetVxStrainRate2d{{{*/ +void TriaInput::GetVxStrainRate2d(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussTria* gauss){ + + /*Intermediary*/ + int numnodes=this->NumberofNodes(); + IssmDouble* B=xNew(numnodes*(NDOF2*numnodes)); + IssmDouble* velocity=xNew(numnodes*NDOF2); + + /*Get B matrix: */ + GetBMacAyeal(&B[0],xyz_list,gauss); + + /*Here, we are computing the strain rate of (vx,0)*/ + for(int i=0;ivalues[i]; + velocity[2*i+1]=0.; + } + /*Get epsilon(vx) = B*velocity*/ + MatrixMultiply(&B[0],numnodes,NDOF2*numnodes,0, + &velocity[0],NDOF2*numnodes,1,0, + epsilonvx,0); + + /*clean up*/ + xDelete(B); + xDelete(velocity); +} +/*}}}*/ +/*FUNCTION TriaInput::GetVyStrainRate2d{{{*/ +void TriaInput::GetVyStrainRate2d(IssmDouble* epsilonvy,IssmDouble* xyz_list, GaussTria* gauss){ + + /*Intermediary*/ + int numnodes=this->NumberofNodes(); + IssmDouble* B=xNew(numnodes*(NDOF2*numnodes)); + IssmDouble* velocity=xNew(numnodes*NDOF2); + + /*Get B matrix: */ + GetBMacAyeal(&B[0],xyz_list,gauss); + + /*Here, we are computing the strain rate of (0,vy)*/ + for(int i=0;ivalues[i]; + } + /*Get epsilon(vy) = B*velocity*/ + MatrixMultiply(&B[0],numnodes,NDOF2*numnodes,0, + &velocity[0],NDOF2*numnodes,1,0, + epsilonvy,0); + + /*clean up*/ + xDelete(B); + xDelete(velocity); +} +/*}}}*/ +/*FUNCTION TriaInput::ChangeEnum{{{*/ +void TriaInput::ChangeEnum(int newenumtype){ + this->enum_type=newenumtype; +} +/*}}}*/ +/*FUNCTION TriaInput::GetInputAverage{{{*/ +void TriaInput::GetInputAverage(IssmDouble* pvalue){ + + int numnodes = this->NumberofNodes(); + IssmDouble numnodesd = reCast(numnodes); + IssmDouble value = 0.; + + for(int i=0;i(1); + times=xNew(1); + + this->GetInputAverage(&outvalues[0]); + times[0]=0.; /*we don't have a time*/ + + *pvalues=outvalues; + *ptimes=times; + *pnumtimes=numtimes; +} +/*}}}*/ +/*FUNCTION TriaInput::GetInputUpToCurrentTimeAverages{{{*/ +void TriaInput::GetInputUpToCurrentTimeAverages(IssmDouble** pvalues, IssmDouble** ptimes, int* pnumtimes, IssmDouble currenttime){ + + IssmDouble* outvalues=NULL; + IssmDouble* times=NULL; + int numtimes; + + /*this is not a transient forcing, so we only have 1 value, steady state: */ + numtimes=1; + outvalues=xNew(1); + times=xNew(1); + + this->GetInputAverage(&outvalues[0]); + times[0]=currenttime; /*we don't have a time*/ + + *pvalues=outvalues; + *ptimes=times; + *pnumtimes=numtimes; +} +/*}}}*/ + +/*Intermediary*/ +/*FUNCTION TriaInput::SquareMin{{{*/ +void TriaInput::SquareMin(IssmDouble* psquaremin,Parameters* parameters){ + + int numnodes=this->NumberofNodes(); + IssmDouble squaremin; + + /*Now, figure out minimum of valuescopy: */ + squaremin=pow(this->values[0],2); + for(int i=1;ivalues[i],2)values[i],2); + } + /*Assign output pointers:*/ + *psquaremin=squaremin; +} +/*}}}*/ +/*FUNCTION TriaInput::ContrainMin{{{*/ +void TriaInput::ConstrainMin(IssmDouble minimum){ + + const int numnodes = this->NumberofNodes(); + for(int i=0;iNumberofNodes(); + + for(int i=0;inorm) norm=fabs(values[i]); + return norm; +} +/*}}}*/ +/*FUNCTION TriaInput::Max{{{*/ +IssmDouble TriaInput::Max(void){ + + const int numnodes=this->NumberofNodes(); + IssmDouble max=values[0]; + + for(int i=1;imax) max=values[i]; + } + return max; +} +/*}}}*/ +/*FUNCTION TriaInput::MaxAbs{{{*/ +IssmDouble TriaInput::MaxAbs(void){ + + const int numnodes=this->NumberofNodes(); + IssmDouble max=fabs(values[0]); + + for(int i=1;imax) max=fabs(values[i]); + } + return max; +} +/*}}}*/ +/*FUNCTION TriaInput::Min{{{*/ +IssmDouble TriaInput::Min(void){ + + const int numnodes=this->NumberofNodes(); + IssmDouble min=values[0]; + + for(int i=1;iNumberofNodes(); + IssmDouble min=fabs(values[0]); + + for(int i=1;iNumberofNodes(); + + for(int i=0;iNumberofNodes(); + IssmDouble noise; + + /*Compute random number between bounds: + * rand() outputs an integer in [0 RAND_MAX] + * (IssmDouble)rand()/RAND_MAX is in [0 1] + */ + noise=min+(max-min)*(IssmDouble)rand()/RAND_MAX; + + for(i=0;iNumberofNodes(); + TriaInput* xtriavertexinput=NULL; + + /*xinput is of the same type, so cast it: */ + if(xinput->ObjectEnum()!=TriaInputEnum) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum())); + xtriavertexinput=(TriaInput*)xinput; + if(xtriavertexinput->element_type!=this->element_type) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum())); + + /*Carry out the AXPY operation depending on type:*/ + for(i=0;ivalues[i]=this->values[i]+scalar*xtriavertexinput->values[i]; + +} +/*}}}*/ +/*FUNCTION TriaInput::Constrain{{{*/ +void TriaInput::Constrain(IssmDouble cm_min, IssmDouble cm_max){ + + int i; + const int numnodes=this->NumberofNodes(); + + if(!xIsNan(cm_min)) for(i=0;ivalues[i]values[i]=cm_min; + if(!xIsNan(cm_max)) for(i=0;ivalues[i]>cm_max)this->values[i]=cm_max; + +} +/*}}}*/ +/*FUNCTION TriaInput::GetVectorFromInputs{{{*/ +void TriaInput::GetVectorFromInputs(Vector* vector,int* doflist){ + + const int numnodes=this->NumberofNodes(); + vector->SetValues(numnodes,doflist,this->values,INS_VAL); + +} /*}}}*/ +/*FUNCTION TriaInput::PointwiseMin{{{*/ +Input* TriaInput::PointwiseMin(Input* inputB){ + + /*Ouput*/ + TriaInput* outinput=NULL; + + /*Intermediaries*/ + int i; + TriaInput *xinputB = NULL; + const int numnodes = this->NumberofNodes(); + IssmDouble *minvalues = xNew(numnodes); + + /*Check that inputB is of the same type*/ + if(inputB->ObjectEnum()!=TriaInputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum())); + xinputB=(TriaInput*)inputB; + if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum())); + + /*Create point wise min*/ + for(i=0;ivalues[i] > xinputB->values[i]) minvalues[i]=xinputB->values[i]; + else minvalues[i]=this->values[i]; + } + + /*Create new Tria vertex input (copy of current input)*/ + outinput=new TriaInput(this->enum_type,&minvalues[0],this->element_type); + + /*Return output pointer*/ + xDelete(minvalues); + return outinput; + +} +/*}}}*/ +/*FUNCTION TriaInput::PointwiseMax{{{*/ +Input* TriaInput::PointwiseMax(Input* inputB){ + + /*Ouput*/ + TriaInput* outinput=NULL; + + /*Intermediaries*/ + int i; + TriaInput *xinputB = NULL; + const int numnodes = this->NumberofNodes(); + IssmDouble *maxvalues = xNew(numnodes); + + /*Check that inputB is of the same type*/ + if(inputB->ObjectEnum()!=TriaInputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum())); + xinputB=(TriaInput*)inputB; + if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum())); + + /*Create point wise max*/ + for(i=0;ivalues[i] < xinputB->values[i]) maxvalues[i]=xinputB->values[i]; + else maxvalues[i]=this->values[i]; + } + + /*Create new Tria vertex input (copy of current input)*/ + outinput=new TriaInput(this->enum_type,&maxvalues[0],this->element_type); + + /*Return output pointer*/ + xDelete(maxvalues); + return outinput; + +} +/*}}}*/ +/*FUNCTION TriaInput::Configure{{{*/ +void TriaInput::Configure(Parameters* parameters){ + /*do nothing: */ +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h (revision 15298) +++ /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h (revision 15298) @@ -0,0 +1,85 @@ +/*! \file TriaInput.h + * \brief: header file for TriaInput object + */ + +#ifndef _TRIAINPUT_H_ +#define _TRIAINPUT_H_ + +/*Headers:*/ +/*{{{*/ +#include "./Input.h" +#include "../Elements/TriaRef.h" +class GaussTria; +class GaussPenta; +/*}}}*/ + +class TriaInput: public Input,public TriaRef{ + + public: + int enum_type; + IssmDouble* values; + + /*TriaInput constructors, destructors*/ + TriaInput(); + TriaInput(int enum_type,IssmDouble* values,int element_type_in=P1Enum); + ~TriaInput(); + + /*Object virtual functions definitions*/ + void Echo(); + void DeepEcho(); + int Id(); + int ObjectEnum(); + Object *copy(); + + /*TriaInput management:*/ + int InstanceEnum(); + Input* SpawnTriaInput(int* indices); + Input* PointwiseDivide(Input* inputB){_error_("not implemented yet");}; + Input* PointwiseMin(Input* inputB); + Input* PointwiseMax(Input* inputB); + ElementResult* SpawnResult(int step, IssmDouble time); + void AddTimeValues(IssmDouble* values,int step,IssmDouble time){_error_("not supported yet");}; + void Configure(Parameters* parameters); + + /*numerics*/ + void GetInputValue(bool* pvalue){_error_("not implemented yet");} + void GetInputValue(int* pvalue){_error_("not implemented yet");} + void GetInputValue(IssmDouble* pvalue){_error_("not implemented yet");} + void GetInputValue(IssmDouble* pvalue,GaussTria* gauss); + void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){_error_("not implemented yet");}; + void GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time){_error_("not implemented yet");}; + void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time){_error_("not implemented yet");}; + void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index){_error_("not implemented yet");}; + void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,int index){_error_("not implemented yet");}; + void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss); + void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss){_error_("not implemented yet");}; + void GetInputAverage(IssmDouble* pvalue); + void GetInputAllTimeAverages(IssmDouble** pvalues,IssmDouble** ptimes, int* pnumtimes); + void GetInputUpToCurrentTimeAverages(IssmDouble** pvalues, IssmDouble** ptimes, int* pnumtimes, IssmDouble currenttime); + void GetVxStrainRate2d(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussTria* gauss); + void GetVyStrainRate2d(IssmDouble* epsilonvy,IssmDouble* xyz_list, GaussTria* gauss); + void GetVxStrainRate3d(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussPenta* gauss){_error_("not implemented yet");}; + void GetVyStrainRate3d(IssmDouble* epsilonvy,IssmDouble* xyz_list, GaussPenta* gauss){_error_("not implemented yet");}; + void GetVzStrainRate3d(IssmDouble* epsilonvz,IssmDouble* xyz_list, GaussPenta* gauss){_error_("not implemented yet");}; + void GetVxStrainRate3dPattyn(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussPenta* gauss){_error_("not implemented yet");}; + void GetVyStrainRate3dPattyn(IssmDouble* epsilonvy,IssmDouble* xyz_list, GaussPenta* gauss){_error_("not implemented yet");}; + void ChangeEnum(int newenumtype); + + void SquareMin(IssmDouble* psquaremin,Parameters* parameters); + void ConstrainMin(IssmDouble minimum); + void Scale(IssmDouble scale_factor); + void ArtificialNoise(IssmDouble min,IssmDouble max); + void AXPY(Input* xinput,IssmDouble scalar); + void Constrain(IssmDouble cm_min, IssmDouble cm_max); + IssmDouble InfinityNorm(void); + IssmDouble Max(void); + IssmDouble MaxAbs(void); + IssmDouble Min(void); + IssmDouble MinAbs(void); + void Extrude(void){_error_("not supported yet");}; + void VerticallyIntegrate(Input* thickness_input){_error_("not supported yet");}; + void GetVectorFromInputs(Vector* vector,int* doflist); + /*}}}*/ + +}; +#endif /* _TRIAINPUT_H */ Index: sm/trunk-jpl/src/c/classes/Inputs/TriaP1Input.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/TriaP1Input.cpp (revision 15297) +++ (revision ) @@ -1,438 +1,0 @@ -/*!\file TriaP1Input.c - * \brief: implementation of the TriaP1Input object - */ - -#ifdef HAVE_CONFIG_H - #include -#else -#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" -#endif - -#include "../classes.h" -#include "../../shared/shared.h" - -/*TriaP1Input constructors and destructor*/ -/*FUNCTION TriaP1Input::TriaP1Input(){{{*/ -TriaP1Input::TriaP1Input(){ - return; -} -/*}}}*/ -/*FUNCTION TriaP1Input::TriaP1Input(int in_enum_type,IssmDouble* values){{{*/ -TriaP1Input::TriaP1Input(int in_enum_type,IssmDouble* in_values) - :TriaRef(1) -{ - - /*Set TriaRef*/ - this->SetElementType(P1Enum,0); - this->element_type=P1Enum; - - /*Set Enum*/ - enum_type=in_enum_type; - - /*Set values*/ - values[0]=in_values[0]; - values[1]=in_values[1]; - values[2]=in_values[2]; -} -/*}}}*/ -/*FUNCTION TriaP1Input::~TriaP1Input(){{{*/ -TriaP1Input::~TriaP1Input(){ - return; -} -/*}}}*/ - -/*Object virtual functions definitions:*/ -/*FUNCTION TriaP1Input::Echo {{{*/ -void TriaP1Input::Echo(void){ - this->DeepEcho(); -} -/*}}}*/ -/*FUNCTION TriaP1Input::DeepEcho{{{*/ -void TriaP1Input::DeepEcho(void){ - - _printf_("TriaP1Input:\n"); - _printf_(" enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")\n"); - _printf_(" values: [" << this->values[0] << " " << this->values[1] << " " << this->values[2] << "]\n"); -} -/*}}}*/ -/*FUNCTION TriaP1Input::Id{{{*/ -int TriaP1Input::Id(void){ return -1; } -/*}}}*/ -/*FUNCTION TriaP1Input::ObjectEnum{{{*/ -int TriaP1Input::ObjectEnum(void){ - - return TriaP1InputEnum; - -} -/*}}}*/ -/*FUNCTION TriaP1Input::copy{{{*/ -Object* TriaP1Input::copy() { - - return new TriaP1Input(this->enum_type,this->values); - -} -/*}}}*/ - -/*TriaP1Input management*/ -/*FUNCTION TriaP1Input::InstanceEnum{{{*/ -int TriaP1Input::InstanceEnum(void){ - - return this->enum_type; - -} -/*}}}*/ -/*FUNCTION TriaP1Input::SpawnTriaInput{{{*/ -Input* TriaP1Input::SpawnTriaInput(int* indices){ - - /*output*/ - TriaP1Input* outinput=NULL; - - /*Create new Tria input (copy of current input)*/ - outinput=new TriaP1Input(this->enum_type,&this->values[0]); - - /*Assign output*/ - return outinput; - -} -/*}}}*/ -/*FUNCTION TriaP1Input::SpawnResult{{{*/ -ElementResult* TriaP1Input::SpawnResult(int step, IssmDouble time){ - - return new TriaP1ElementResult(this->enum_type,this->values,step,time); - -} -/*}}}*/ - -/*Object functions*/ -/*FUNCTION TriaP1Input::GetInputValue(IssmDouble* pvalue,GaussTria* gauss){{{*/ -void TriaP1Input::GetInputValue(IssmDouble* pvalue,GaussTria* gauss){ - - /*Call TriaRef function*/ - TriaRef::GetInputValue(pvalue,&values[0],gauss); - -} -/*}}}*/ -/*FUNCTION TriaP1Input::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, GaussTria* gauss){{{*/ -void TriaP1Input::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, GaussTria* gauss){ - - /*Call TriaRef function*/ - TriaRef::GetInputDerivativeValue(p,&values[0],xyz_list,gauss); -} -/*}}}*/ -/*FUNCTION TriaP1Input::GetVxStrainRate2d{{{*/ -void TriaP1Input::GetVxStrainRate2d(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussTria* gauss){ - - /*Intermediary*/ - int i; - const int numnodes=3; - IssmDouble B[3][NDOF2*numnodes]; - IssmDouble velocity[3][NDOF2]; - - /*Get B matrix: */ - GetBMacAyeal(&B[0][0], xyz_list, gauss); - - /*Here, we are computing the strain rate of (vx,0)*/ - for(i=0;i<3;i++){ - velocity[i][0]=this->values[i]; - velocity[i][1]=0.0; - } - /*Get epsilon(vx) = B*velocity*/ - MatrixMultiply( &B[0][0],3,NDOF2*numnodes,0, - &velocity[0][0],NDOF2*numnodes,1,0, - epsilonvx,0); -} -/*}}}*/ -/*FUNCTION TriaP1Input::GetVyStrainRate2d{{{*/ -void TriaP1Input::GetVyStrainRate2d(IssmDouble* epsilonvy,IssmDouble* xyz_list, GaussTria* gauss){ - - /*Intermediary*/ - int i; - const int numnodes=3; - IssmDouble B[3][NDOF2*numnodes]; - IssmDouble velocity[3][NDOF2]; - - /*Get B matrix: */ - GetBMacAyeal(&B[0][0], xyz_list, gauss); - - /*Here, we are computing the strain rate of (0,vy)*/ - for(i=0;i<3;i++){ - velocity[i][0]=0.0; - velocity[i][1]=this->values[i]; - } - /*Get epsilon(vy) = B*velocity*/ - MatrixMultiply( &B[0][0],3,NDOF2*numnodes,0, - &velocity[0][0],NDOF2*numnodes,1,0, - epsilonvy,0); -} -/*}}}*/ -/*FUNCTION TriaP1Input::ChangeEnum{{{*/ -void TriaP1Input::ChangeEnum(int newenumtype){ - this->enum_type=newenumtype; -} -/*}}}*/ -/*FUNCTION TriaP1Input::GetInputAverage{{{*/ -void TriaP1Input::GetInputAverage(IssmDouble* pvalue){ - *pvalue=1./3.*(values[0]+values[1]+values[2]); -} -/*}}}*/ -/*FUNCTION TriaP1Input::GetInputAllTimeAverages{{{*/ -void TriaP1Input::GetInputAllTimeAverages(IssmDouble** pvalues,IssmDouble** ptimes, int* pnumtimes){ - - IssmDouble* outvalues=NULL; - IssmDouble* times=NULL; - int numtimes; - - /*this is not a transient forcing, so we only have 1 value, steady state: */ - numtimes=1; - outvalues=xNew(1); - times=xNew(1); - - outvalues[0]=1./3.*(values[0]+values[1]+values[2]); - times[0]=0; /*we don't have a time*/ - - *pvalues=outvalues; - *ptimes=times; - *pnumtimes=numtimes; -} -/*}}}*/ -/*FUNCTION TriaP1Input::GetInputUpToCurrentTimeAverages{{{*/ -void TriaP1Input::GetInputUpToCurrentTimeAverages(IssmDouble** pvalues, IssmDouble** ptimes, int* pnumtimes, IssmDouble currenttime){ - - IssmDouble* outvalues=NULL; - IssmDouble* times=NULL; - int numtimes; - - /*this is not a transient forcing, so we only have 1 value, steady state: */ - numtimes=1; - outvalues=xNew(1); - times=xNew(1); - - outvalues[0]=1./3.*(values[0]+values[1]+values[2]); - times[0]=currenttime; /*we don't have a time*/ - - *pvalues=outvalues; - *ptimes=times; - *pnumtimes=numtimes; -} -/*}}}*/ - -/*Intermediary*/ -/*FUNCTION TriaP1Input::SquareMin{{{*/ -void TriaP1Input::SquareMin(IssmDouble* psquaremin,Parameters* parameters){ - - int i; - const int numnodes=3; - IssmDouble valuescopy[numnodes]; - IssmDouble squaremin; - - /*First, copy values, to process units if requested: */ - for(i=0;ivalues[i]; - - /*Now, figure out minimum of valuescopy: */ - squaremin=pow(valuescopy[0],2); - for(i=1;inorm) norm=fabs(values[i]); - return norm; -} -/*}}}*/ -/*FUNCTION TriaP1Input::Max{{{*/ -IssmDouble TriaP1Input::Max(void){ - - const int numnodes=3; - IssmDouble max=values[0]; - - for(int i=1;imax) max=values[i]; - } - return max; -} -/*}}}*/ -/*FUNCTION TriaP1Input::MaxAbs{{{*/ -IssmDouble TriaP1Input::MaxAbs(void){ - - const int numnodes=3; - IssmDouble max=fabs(values[0]); - - for(int i=1;imax) max=fabs(values[i]); - } - return max; -} -/*}}}*/ -/*FUNCTION TriaP1Input::Min{{{*/ -IssmDouble TriaP1Input::Min(void){ - - const int numnodes=3; - IssmDouble min=values[0]; - - for(int i=1;iObjectEnum()){ - - case TriaP1InputEnum : - for(i=0;ivalues[i]=this->values[i]+scalar*xtriavertexinput->values[i]; - return; - - default : - _error_("not implemented yet"); - } - -} -/*}}}*/ -/*FUNCTION TriaP1Input::Constrain{{{*/ -void TriaP1Input::Constrain(IssmDouble cm_min, IssmDouble cm_max){ - - int i; - const int numnodes=3; - - if(!xIsNan(cm_min)) for(i=0;ivalues[i]values[i]=cm_min; - if(!xIsNan(cm_max)) for(i=0;ivalues[i]>cm_max)this->values[i]=cm_max; - -} -/*}}}*/ -/*FUNCTION TriaP1Input::GetVectorFromInputs{{{*/ -void TriaP1Input::GetVectorFromInputs(Vector* vector,int* doflist){ - - const int numvertices=3; - vector->SetValues(numvertices,doflist,this->values,INS_VAL); - -} /*}}}*/ -/*FUNCTION TriaP1Input::PointwiseMin{{{*/ -Input* TriaP1Input::PointwiseMin(Input* inputB){ - - /*Ouput*/ - TriaP1Input* outinput=NULL; - - /*Intermediaries*/ - int i; - TriaP1Input *xinputB = NULL; - const int numnodes = 3; - IssmDouble minvalues[numnodes]; - - /*Check that inputB is of the same type*/ - if (inputB->ObjectEnum()!=TriaP1InputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum())); - xinputB=(TriaP1Input*)inputB; - - /*Create point wise min*/ - for(i=0;ivalues[i] > xinputB->values[i]) minvalues[i]=xinputB->values[i]; - else minvalues[i]=this->values[i]; - } - - /*Create new Tria vertex input (copy of current input)*/ - outinput=new TriaP1Input(this->enum_type,&minvalues[0]); - - /*Return output pointer*/ - return outinput; - -} -/*}}}*/ -/*FUNCTION TriaP1Input::PointwiseMax{{{*/ -Input* TriaP1Input::PointwiseMax(Input* inputB){ - - /*Ouput*/ - TriaP1Input* outinput=NULL; - - /*Intermediaries*/ - int i; - TriaP1Input *xinputB = NULL; - const int numnodes = 3; - IssmDouble maxvalues[numnodes]; - - /*Check that inputB is of the same type*/ - if (inputB->ObjectEnum()!=TriaP1InputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum())); - xinputB=(TriaP1Input*)inputB; - - /*Create point wise max*/ - for(i=0;ivalues[i] < xinputB->values[i]) maxvalues[i]=xinputB->values[i]; - else maxvalues[i]=this->values[i]; - } - - /*Create new Tria vertex input (copy of current input)*/ - outinput=new TriaP1Input(this->enum_type,&maxvalues[0]); - - /*Return output pointer*/ - return outinput; - -} -/*}}}*/ -/*FUNCTION TriaP1Input::Configure{{{*/ -void TriaP1Input::Configure(Parameters* parameters){ - /*do nothing: */ -} -/*}}}*/ Index: sm/trunk-jpl/src/c/classes/Inputs/TriaP1Input.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/TriaP1Input.h (revision 15297) +++ (revision ) @@ -1,86 +1,0 @@ -/*! \file TriaP1Input.h - * \brief: header file for TriaP1Input object - */ - -#ifndef _TRIAP1INPUT_H_ -#define _TRIAP1INPUT_H_ - -/*Headers:*/ -/*{{{*/ -#include "./Input.h" -#include "../Elements/TriaRef.h" -class GaussTria; -class GaussPenta; -/*}}}*/ - -class TriaP1Input: public Input,public TriaRef{ - - public: - /*just hold 3 values for 3 vertices: */ - int enum_type; - IssmDouble values[3]; - - /*TriaP1Input constructors, destructors: {{{*/ - TriaP1Input(); - TriaP1Input(int enum_type,IssmDouble* values); - ~TriaP1Input(); - /*}}}*/ - /*Object virtual functions definitions:{{{ */ - void Echo(); - void DeepEcho(); - int Id(); - int ObjectEnum(); - Object *copy(); - /*}}}*/ - /*TriaP1Input management: {{{*/ - int InstanceEnum(); - Input* SpawnTriaInput(int* indices); - Input* PointwiseDivide(Input* inputB){_error_("not implemented yet");}; - Input* PointwiseMin(Input* inputB); - Input* PointwiseMax(Input* inputB); - ElementResult* SpawnResult(int step, IssmDouble time); - void AddTimeValues(IssmDouble* values,int step,IssmDouble time){_error_("not supported yet");}; - void Configure(Parameters* parameters); - /*}}}*/ - /*numerics: {{{*/ - void GetInputValue(bool* pvalue){_error_("not implemented yet");} - void GetInputValue(int* pvalue){_error_("not implemented yet");} - void GetInputValue(IssmDouble* pvalue){_error_("not implemented yet");} - void GetInputValue(IssmDouble* pvalue,GaussTria* gauss); - void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){_error_("not implemented yet");}; - void GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time){_error_("not implemented yet");}; - void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time){_error_("not implemented yet");}; - void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index){_error_("not implemented yet");}; - void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,int index){_error_("not implemented yet");}; - void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss); - void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss){_error_("not implemented yet");}; - void GetInputAverage(IssmDouble* pvalue); - void GetInputAllTimeAverages(IssmDouble** pvalues,IssmDouble** ptimes, int* pnumtimes); - void GetInputUpToCurrentTimeAverages(IssmDouble** pvalues, IssmDouble** ptimes, int* pnumtimes, IssmDouble currenttime); - void GetVxStrainRate2d(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussTria* gauss); - void GetVyStrainRate2d(IssmDouble* epsilonvy,IssmDouble* xyz_list, GaussTria* gauss); - void GetVxStrainRate3d(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussPenta* gauss){_error_("not implemented yet");}; - void GetVyStrainRate3d(IssmDouble* epsilonvy,IssmDouble* xyz_list, GaussPenta* gauss){_error_("not implemented yet");}; - void GetVzStrainRate3d(IssmDouble* epsilonvz,IssmDouble* xyz_list, GaussPenta* gauss){_error_("not implemented yet");}; - void GetVxStrainRate3dPattyn(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussPenta* gauss){_error_("not implemented yet");}; - void GetVyStrainRate3dPattyn(IssmDouble* epsilonvy,IssmDouble* xyz_list, GaussPenta* gauss){_error_("not implemented yet");}; - void ChangeEnum(int newenumtype); - - void SquareMin(IssmDouble* psquaremin,Parameters* parameters); - void ConstrainMin(IssmDouble minimum); - void Scale(IssmDouble scale_factor); - void ArtificialNoise(IssmDouble min,IssmDouble max); - void AXPY(Input* xinput,IssmDouble scalar); - void Constrain(IssmDouble cm_min, IssmDouble cm_max); - IssmDouble InfinityNorm(void); - IssmDouble Max(void); - IssmDouble MaxAbs(void); - IssmDouble Min(void); - IssmDouble MinAbs(void); - void Extrude(void){_error_("not supported yet");}; - void VerticallyIntegrate(Input* thickness_input){_error_("not supported yet");}; - void GetVectorFromInputs(Vector* vector,int* doflist); - /*}}}*/ - -}; -#endif /* _TRIAP1INPUT_H */ Index: /issm/trunk-jpl/src/c/classes/Loads/Icefront.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Icefront.cpp (revision 15297) +++ /issm/trunk-jpl/src/c/classes/Loads/Icefront.cpp (revision 15298) @@ -41,5 +41,6 @@ int segment_width; int element; - int num_nodes; + int numnodes; + int numvertices; int dim; int numberofelements; @@ -88,8 +89,12 @@ else _error_("in_icefront_type " << EnumToStringx(in_icefront_type) << " not supported yet!"); - if (in_icefront_type==PattynIceFrontEnum || in_icefront_type==StokesIceFrontEnum) - num_nodes=4; - else - num_nodes=2; + if (in_icefront_type==PattynIceFrontEnum || in_icefront_type==StokesIceFrontEnum){ + numnodes=4; + numvertices=4; + } + else{ + numnodes=2; + numvertices=2; + } /*Fill*/ @@ -101,6 +106,6 @@ /*Hooks: */ - this->hnodes=new Hook(icefront_node_ids,num_nodes); - this->hvertices=new Hook(icefront_vertex_ids,num_nodes); + this->hnodes=new Hook(icefront_node_ids,numnodes); + this->hvertices=new Hook(icefront_vertex_ids,numvertices); this->helement=new Hook(&icefront_eid,1); this->hmatpar=new Hook(&icefront_mparid,1); @@ -246,9 +251,7 @@ /*Checks in debugging mode*/ - /*{{{*/ _assert_(nodes); _assert_(element); _assert_(matpar); - /*}}}*/ /*Retrieve parameters: */ @@ -496,23 +499,25 @@ ElementVector* Icefront::CreatePVectorDiagnosticMacAyeal2d(void){ - /*Constants*/ - const int numnodes= NUMVERTICESSEG; - const int numdofs = numnodes *NDOF2; - /*Intermediary*/ - int ig,index1,index2,fill; - IssmDouble Jdet; - IssmDouble thickness,bed,pressure,ice_pressure,rho_water,rho_ice,gravity; - IssmDouble water_pressure,air_pressure,surface_under_water,base_under_water; - IssmDouble xyz_list[numnodes][3]; - IssmDouble normal[2]; - IssmDouble L[2]; + int ig,index1,index2,fill; + IssmDouble Jdet; + IssmDouble thickness,bed,pressure,ice_pressure,rho_water,rho_ice,gravity; + IssmDouble water_pressure,air_pressure,surface_under_water,base_under_water; + IssmDouble xyz_list[NUMVERTICESSEG][3]; + IssmDouble normal[2]; GaussTria *gauss; + /*return of element is on water*/ Tria* tria=((Tria*)element); - - /*Initialize Element vector and return if necessary*/ if(tria->IsOnWater()) return NULL; + + /*Fetch number of nodes and dof for this finite element*/ + //int numnodes = this->NumberofNodes(); + int numnodes = 2; + int numdof = numnodes*NDOF2; + + /*Initialize Element vector and vectors*/ ElementVector* pe=new ElementVector(nodes,NUMVERTICESSEG,this->parameters,MacAyealApproximationEnum); + IssmDouble* basis = xNew(numnodes); /*Retrieve all inputs and parameters*/ @@ -558,9 +563,9 @@ tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss); - tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2); + tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2); for (int i=0;ivalues[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*L[i]; - pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*L[i]; + pe->values[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i]; + pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i]; } } @@ -570,4 +575,5 @@ /*Clean up and return*/ + xDelete(basis); delete gauss; return pe; Index: /issm/trunk-jpl/src/c/classes/Materials/Matdamageice.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matdamageice.cpp (revision 15297) +++ /issm/trunk-jpl/src/c/classes/Materials/Matdamageice.cpp (revision 15298) @@ -630,5 +630,5 @@ IssmDouble values[3]; for (int i=0;i<3;i++) values[i]=vector[((Tria*)element)->vertices[i]->Pid()]; - this->inputs->AddInput(new TriaP1Input(name,values)); + this->inputs->AddInput(new TriaInput(name,values)); return; } @@ -672,5 +672,5 @@ IssmDouble values[3]; for (int i=0;i<3;i++) values[i]=vector[((Tria*)element)->vertices[i]->Sid()]; //use sid list, to index into serial oriented vector - this->inputs->AddInput(new TriaP1Input(name,values)); + this->inputs->AddInput(new TriaInput(name,values)); /*Special case for rheology B in 2D: Pourave land for this solution{{{*/ if(name==MaterialsRheologyBEnum){ @@ -685,5 +685,5 @@ if(dim==2){ /*Dupliacte rheology input: */ - this->inputs->AddInput(new TriaP1Input(MaterialsRheologyBbarEnum,values)); + this->inputs->AddInput(new TriaInput(MaterialsRheologyBbarEnum,values)); } } @@ -759,5 +759,5 @@ if (iomodel->Data(MaterialsRheologyBEnum)) { for(i=0;iData(MaterialsRheologyBEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)]; - this->inputs->AddInput(new TriaP1Input(MaterialsRheologyBbarEnum,nodeinputs)); + this->inputs->AddInput(new TriaInput(MaterialsRheologyBbarEnum,nodeinputs)); } @@ -765,5 +765,5 @@ if (iomodel->Data(MaterialsRheologyNEnum)) { for(i=0;iData(MaterialsRheologyNEnum)[index]; - this->inputs->AddInput(new TriaP1Input(MaterialsRheologyNEnum,nodeinputs)); + this->inputs->AddInput(new TriaInput(MaterialsRheologyNEnum,nodeinputs)); } @@ -771,5 +771,5 @@ if (iomodel->Data(MaterialsRheologyZEnum)) { for(i=0;iData(MaterialsRheologyZEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)]; - this->inputs->AddInput(new TriaP1Input(MaterialsRheologyZbarEnum,nodeinputs)); + this->inputs->AddInput(new TriaInput(MaterialsRheologyZbarEnum,nodeinputs)); } @@ -785,5 +785,5 @@ for(j=0;jData(InversionMinParametersEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; for(j=0;jData(InversionMaxParametersEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; - this->inputs->AddInput(new ControlInput(MaterialsRheologyBbarEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(MaterialsRheologyBbarEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); } break; @@ -794,5 +794,5 @@ for(j=0;jData(InversionMinParametersEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; for(j=0;jData(InversionMaxParametersEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; - this->inputs->AddInput(new ControlInput(MaterialsRheologyZbarEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(MaterialsRheologyZbarEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); } break; Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 15297) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 15298) @@ -13,5 +13,5 @@ #include "../Inputs/Input.h" #include "../Inputs/Inputs.h" -#include "../Inputs/TriaP1Input.h" +#include "../Inputs/TriaInput.h" #include "../Inputs/PentaP1Input.h" #include "../Inputs/ControlInput.h" @@ -564,5 +564,5 @@ IssmDouble values[3]; for (int i=0;i<3;i++) values[i]=vector[((Tria*)element)->vertices[i]->Pid()]; - this->inputs->AddInput(new TriaP1Input(name,values)); + this->inputs->AddInput(new TriaInput(name,values)); return; } @@ -606,5 +606,5 @@ IssmDouble values[3]; for (int i=0;i<3;i++) values[i]=vector[((Tria*)element)->vertices[i]->Sid()]; //index into serial oriented vector - this->inputs->AddInput(new TriaP1Input(name,values)); + this->inputs->AddInput(new TriaInput(name,values)); /*Special case for rheology B in 2D: Pourave land for this solution{{{*/ if(name==MaterialsRheologyBEnum){ @@ -619,5 +619,5 @@ if(dim==2){ /*Dupliacte rheology input: */ - this->inputs->AddInput(new TriaP1Input(MaterialsRheologyBbarEnum,values)); + this->inputs->AddInput(new TriaInput(MaterialsRheologyBbarEnum,values)); } } @@ -699,5 +699,5 @@ if (iomodel->Data(MaterialsRheologyBEnum)) { for(i=0;iData(MaterialsRheologyBEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)]; - this->inputs->AddInput(new TriaP1Input(MaterialsRheologyBbarEnum,nodeinputs)); + this->inputs->AddInput(new TriaInput(MaterialsRheologyBbarEnum,nodeinputs)); } @@ -705,5 +705,5 @@ if (iomodel->Data(MaterialsRheologyNEnum)) { for(i=0;iData(MaterialsRheologyNEnum)[index]; - this->inputs->AddInput(new TriaP1Input(MaterialsRheologyNEnum,nodeinputs)); + this->inputs->AddInput(new TriaInput(MaterialsRheologyNEnum,nodeinputs)); } @@ -719,5 +719,5 @@ for(j=0;jData(InversionMinParametersEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; for(j=0;jData(InversionMaxParametersEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; - this->inputs->AddInput(new ControlInput(MaterialsRheologyBbarEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + this->inputs->AddInput(new ControlInput(MaterialsRheologyBbarEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); } break; Index: /issm/trunk-jpl/src/c/classes/classes.h =================================================================== --- /issm/trunk-jpl/src/c/classes/classes.h (revision 15297) +++ /issm/trunk-jpl/src/c/classes/classes.h (revision 15298) @@ -58,5 +58,5 @@ #include "./Inputs/IntInput.h" #include "./Inputs/PentaP1Input.h" -#include "./Inputs/TriaP1Input.h" +#include "./Inputs/TriaInput.h" #include "./Inputs/ControlInput.h" #include "./Inputs/DatasetInput.h" Index: /issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp (revision 15297) +++ /issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp (revision 15298) @@ -388,16 +388,40 @@ /*update static arrays*/ switch(iv){ - case 0: - coord1=1; coord2=0; coord3=0; + case 0: coord1=1.; coord2=0.; coord3=0.; break; + case 1: coord1=0.; coord2=1.; coord3=0.; break; + case 2: coord1=0.; coord2=0.; coord3=1.; break; + default: _error_("vertex index should be in [0 2]"); + } + +} +/*}}}*/ +/*FUNCTION GaussTria::GaussNode{{{*/ +void GaussTria::GaussNode(int numnodes,int iv){ + + /*in debugging mode: check that the default constructor has been called*/ + _assert_(numgauss==-1); + + /*update static arrays*/ + switch(numnodes){ + case 3: //P1 Lagrange element + switch(iv){ + case 0: coord1=1.; coord2=0.; coord3=0.; break; + case 1: coord1=0.; coord2=1.; coord3=0.; break; + case 2: coord1=0.; coord2=0.; coord3=1.; break; + default: _error_("node index should be in [0 2]"); + } break; - case 1: - coord1=0; coord2=1; coord3=0; + case 6: //P2 Lagrange element + switch(iv){ + case 0: coord1=1.; coord2=0.; coord3=0.; break; + case 1: coord1=0.; coord2=1.; coord3=0.; break; + case 2: coord1=0.; coord2=0.; coord3=1.; break; + case 3: coord1=0.; coord2=.5; coord3=.5; break; + case 4: coord1=.5; coord2=0.; coord3=.5; break; + case 5: coord1=.5; coord2=.5; coord3=0.; break; + default: _error_("node index should be in [0 5]"); + } break; - case 2: - coord1=0; coord2=0; coord3=1; - break; - default: - _error_("vertex index should be in [0 2]"); - + default: _error_("supported number of nodes are 3 and 6"); } Index: /issm/trunk-jpl/src/c/classes/gauss/GaussTria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/GaussTria.h (revision 15297) +++ /issm/trunk-jpl/src/c/classes/gauss/GaussTria.h (revision 15298) @@ -40,4 +40,5 @@ void GaussPoint(int ig); void GaussVertex(int iv); + void GaussNode(int numnodes,int iv); void GaussCenter(void); void GaussEdgeCenter(int index1,int index2); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp (revision 15297) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp (revision 15298) @@ -51,4 +51,31 @@ } + if(false){ + /*Quadratic element*/ + int numberofedges; + int element1,element2; + bool my_edge; + + iomodel->Constant(&numberofedges,MeshNumberofedgesEnum); + + for (i=0;i(iomodel->Data(MeshEdgesEnum)[4*i+2])-1; //edges are [node1 node2 elem1 elem2] + element2=reCast(iomodel->Data(MeshEdgesEnum)[4*i+3])-1; //edges are [node1 node2 elem1 elem2] + + /*Check whether we should include this edge (element2 is -2 for boundary edges)*/ + my_edge = iomodel->my_elements[element1]; + if(!my_edge && element2>=0){ + my_edge = iomodel->my_elements[element2]; + } + + /*Add node on edge*/ + if(my_edge){ + nodes->AddObject(new Node(iomodel->nodecounter+numberofvertices+i+1,numberofvertices+i,numberofvertices+i+1,i,iomodel,DiagnosticHorizAnalysisEnum)); + } + } + } + /*Clean fetched data: */ iomodel->DeleteData(9,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBordermacayealEnum,FlowequationBorderstokesEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 15297) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 15298) @@ -384,5 +384,5 @@ StringParamEnum, TriaEnum, - TriaP1InputEnum, + TriaInputEnum, VertexEnum, VertexPIdEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 15297) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 15298) @@ -384,5 +384,5 @@ case StringParamEnum : return "StringParam"; case TriaEnum : return "Tria"; - case TriaP1InputEnum : return "TriaP1Input"; + case TriaInputEnum : return "TriaInput"; case VertexEnum : return "Vertex"; case VertexPIdEnum : return "VertexPId"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 15297) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 15298) @@ -393,5 +393,5 @@ else if (strcmp(name,"StringParam")==0) return StringParamEnum; else if (strcmp(name,"Tria")==0) return TriaEnum; - else if (strcmp(name,"TriaP1Input")==0) return TriaP1InputEnum; + else if (strcmp(name,"TriaInput")==0) return TriaInputEnum; else if (strcmp(name,"Vertex")==0) return VertexEnum; else if (strcmp(name,"VertexPId")==0) return VertexPIdEnum; Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 15297) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 15298) @@ -5147,17 +5147,17 @@ return StringToEnum('Tria')[0] -def TriaP1InputEnum(): - """ - TRIAP1INPUTENUM - Enum of TriaP1Input - - WARNING: DO NOT MODIFY THIS FILE - this file has been automatically generated by src/c/shared/Enum/Synchronize.sh - Please read src/c/shared/Enum/README for more information - - Usage: - macro=TriaP1InputEnum() - """ - - return StringToEnum('TriaP1Input')[0] +def TriaInputEnum(): + """ + TRIAINPUTENUM - Enum of TriaInput + + WARNING: DO NOT MODIFY THIS FILE + this file has been automatically generated by src/c/shared/Enum/Synchronize.sh + Please read src/c/shared/Enum/README for more information + + Usage: + macro=TriaInputEnum() + """ + + return StringToEnum('TriaInput')[0] def VertexEnum(): Index: /issm/trunk-jpl/src/m/enum/TriaInputEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/TriaInputEnum.m (revision 15298) +++ /issm/trunk-jpl/src/m/enum/TriaInputEnum.m (revision 15298) @@ -0,0 +1,11 @@ +function macro=TriaInputEnum() +%TRIAINPUTENUM - Enum of TriaInput +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=TriaInputEnum() + +macro=StringToEnum('TriaInput');