Index: /issm/trunk/src/c/objects/Loads/Icefront.cpp =================================================================== --- /issm/trunk/src/c/objects/Loads/Icefront.cpp (revision 5910) +++ /issm/trunk/src/c/objects/Loads/Icefront.cpp (revision 5911) @@ -318,7 +318,9 @@ /*Checks in debugging mode*/ + /*{{{2*/ ISSMASSERT(nodes); ISSMASSERT(element); ISSMASSERT(matpar); + /*}}}*/ /*Retrieve parameters: */ Index: /issm/trunk/src/c/objects/Loads/Numericalflux.cpp =================================================================== --- /issm/trunk/src/c/objects/Loads/Numericalflux.cpp (revision 5910) +++ /issm/trunk/src/c/objects/Loads/Numericalflux.cpp (revision 5911) @@ -332,4 +332,5 @@ int type; int analysis_type; + ElementMatrix* Ke=NULL; /*recover some parameters*/ @@ -337,11 +338,26 @@ this->parameters->FindParam(&analysis_type,AnalysisTypeEnum); - if (analysis_type==PrognosticAnalysisEnum || analysis_type==BalancedthicknessAnalysisEnum || analysis_type==AdjointBalancedthicknessAnalysisEnum){ - if (type==InternalEnum) CreateKMatrixInternal(Kgg); - else if (type==BoundaryEnum) CreateKMatrixBoundary(Kgg); - else ISSMERROR("type not supported yet"); - } - else{ - ISSMERROR("analysis_type %s not supported yet",EnumToString(analysis_type)); + /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ + switch(analysis_type){ + case PrognosticAnalysisEnum: case BalancedthicknessAnalysisEnum: case AdjointBalancedthicknessAnalysisEnum: + switch(type){ + case InternalEnum: + Ke=CreateKMatrixInternal(); + break; + case BoundaryEnum: + Ke=CreateKMatrixBoundary(); + break; + default: + ISSMERROR("type not supported yet"); + } + break; + default: + ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type)); + } + + /*Add to global matrix*/ + if(Ke){ + Ke->AddToGlobal(Kgg,Kff,Kfs); + delete Ke; } @@ -353,4 +369,5 @@ int type; int analysis_type; + ElementVector* pe=NULL; /*recover some parameters*/ @@ -358,16 +375,26 @@ this->parameters->FindParam(&analysis_type,AnalysisTypeEnum); - if (analysis_type==PrognosticAnalysisEnum || analysis_type==BalancedthicknessAnalysisEnum){ - if (type==InternalEnum) CreatePVectorInternal(pg); - else if (type==BoundaryEnum) CreatePVectorBoundary(pg); - else ISSMERROR("type not supported yet"); - } - else if (analysis_type==AdjointBalancedthicknessAnalysisEnum){ - return; - } - else{ - ISSMERROR("analysis_type %s not supported yet",EnumToString(analysis_type)); - } - + switch(analysis_type){ + case PrognosticAnalysisEnum: case BalancedthicknessAnalysisEnum: case AdjointBalancedthicknessAnalysisEnum: + switch(type){ + case InternalEnum: + pe=CreatePVectorInternal(); + break; + case BoundaryEnum: + pe=CreatePVectorBoundary(); + break; + default: + ISSMERROR("type not supported yet"); + } + break; + default: + ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type)); + } + + /*Add to global matrix*/ + if(pe){ + pe->AddToGlobal(pg,pf); + delete pe; + } } @@ -398,5 +425,5 @@ /*Numericalflux management*/ /*FUNCTION Numericalflux::CreateKMatrixInternal {{{1*/ -void Numericalflux::CreateKMatrixInternal(Mat Kgg){ +ElementMatrix* Numericalflux::CreateKMatrixInternal(void){ /* constants*/ @@ -412,13 +439,13 @@ double Ke_g1[numdof][numdof]; double Ke_g2[numdof][numdof]; - double Ke[numdof][numdof] = {0.0}; - int *doflist = NULL; GaussTria *gauss; + /*Initialize Element matrix and return if necessary*/ + Tria* tria=(Tria*)element; + if(tria->IsOnWater()) return NULL; + ElementMatrix* Ke=NewElementMatrix(nodes,NUMVERTICES_INTERNAL,this->parameters); + + /*Retrieve all inputs and parameters*/ GetVerticesCoordinates(&xyz_list[0][0], nodes,NUMVERTICES_INTERNAL); - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - - /*Retrieve all inputs and parameters we will be needing: */ - Tria* tria=(Tria*)element; Input* vxaverage_input=tria->inputs->GetInput(VxEnum); Input* vyaverage_input=tria->inputs->GetInput(VyEnum); @@ -461,18 +488,15 @@ &Ke_g2[0][0],0); - for(i=0;ivalues[i*numdof+j]+=Ke_g1[i][j]; + for(i=0;ivalues[i*numdof+j]+=Ke_g2[i][j]; + } - /*Free ressources:*/ + /*Clean up and return*/ delete gauss; - xfree((void**)&doflist); - + return Ke; } /*}}}*/ /*FUNCTION Numericalflux::CreateKMatrixBoundary {{{1*/ -void Numericalflux::CreateKMatrixBoundary(Mat Kgg){ +ElementMatrix* Numericalflux::CreateKMatrixBoundary(void){ /* constants*/ @@ -486,12 +510,13 @@ double L[numdof]; double Ke_g[numdof][numdof]; - double Ke[numdof][numdof] = {0.0}; - int *doflist = NULL; GaussTria *gauss; + /*Initialize Element matrix and return if necessary*/ + Tria* tria=(Tria*)element; + if(tria->IsOnWater()) return NULL; + ElementMatrix* Ke=NewElementMatrix(nodes,NUMVERTICES_BOUNDARY,this->parameters); + + /*Retrieve all inputs and parameters*/ GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES_BOUNDARY); - - /*Retrieve all inputs and parameters we will be needing: */ - Tria* tria=(Tria*)element; Input* vxaverage_input=tria->inputs->GetInput(VxEnum); Input* vyaverage_input=tria->inputs->GetInput(VyEnum); @@ -520,8 +545,6 @@ if (UdotN<=0){ /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/ - return; - } - - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + return NULL; + } /* Start looping on the number of gaussian points: */ @@ -544,24 +567,22 @@ &Ke_g[0][0],0); - for(i=0;ivalues[i*numdof+j]+=Ke_g[i][j]; } - MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES); - - /*Free ressources:*/ + /*Clean up and return*/ delete gauss; - xfree((void**)&doflist); + return Ke; } /*}}}*/ /*FUNCTION Numericalflux::CreatePVectorInternal{{{1*/ -void Numericalflux::CreatePVectorInternal(Vec pg){ +ElementVector* Numericalflux::CreatePVectorInternal(void){ /*Nothing added to PVector*/ - return; + return NULL; } /*}}}*/ /*FUNCTION Numericalflux::CreatePVectorBoundary{{{1*/ -void Numericalflux::CreatePVectorBoundary(Vec pg){ +ElementVector* Numericalflux::CreatePVectorBoundary(void){ /* constants*/ @@ -574,12 +595,13 @@ double normal[2]; double L[numdof]; - double pe[numdof] = {0.0}; - int *doflist = NULL; GaussTria *gauss; + /*Initialize Load Vectorand return if necessary*/ + Tria* tria=(Tria*)element; + if(tria->IsOnWater()) return NULL; + ElementVector* pe=NewElementVector(nodes,NUMVERTICES_BOUNDARY,this->parameters); + + /*Retrieve all inputs and parameters*/ GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES_BOUNDARY); - - /*Retrieve all inputs and parameters we will be needing: */ - Tria* tria=(Tria*)element; Input* vxaverage_input=tria->inputs->GetInput(VxEnum); ISSMASSERT(vxaverage_input); Input* vyaverage_input=tria->inputs->GetInput(VyEnum); ISSMASSERT(vyaverage_input); @@ -614,8 +636,6 @@ if (UdotN>0){ /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/ - return; - } - - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + return NULL; + } /* Start looping on the number of gaussian points: */ @@ -634,56 +654,10 @@ DL= - gauss->weight*Jdet*dt*UdotN*thickness; - for(i=0;ivalues[i] += DL*L[i]; + } + + /*Clean up and return*/ delete gauss; - xfree((void**)&doflist); -} -/*}}}*/ -/*FUNCTION Numericalflux::GetDofList {{{1*/ -void Numericalflux::GetDofList(int** pdoflist,int approximation,int setenum){ - - int i,j; - int numberofdofs=0; - int count=0; - int type; - int numberofnodes; - int* doflist=NULL; - - /*Some checks for debugging*/ - ISSMASSERT(nodes); - - /*How many nodes? :*/ - inputs->GetParameterValue(&type,TypeEnum); - switch(type){ - case InternalEnum: - numberofnodes=NUMVERTICES_INTERNAL; break; - case BoundaryEnum: - numberofnodes=NUMVERTICES_BOUNDARY; break; - default: - ISSMERROR("type %s not supported yet",type); - } - - /*Figure out size of doflist: */ - for(i=0;iGetNumberOfDofs(approximation,setenum); - } - - /*Allocate: */ - doflist=(int*)xmalloc(numberofdofs*sizeof(int)); - - /*Populate: */ - count=0; - for(i=0;iGetDofList(doflist+count,approximation,setenum); - count+=nodes[i]->GetNumberOfDofs(approximation,setenum); - } - - /*Assign output pointers:*/ - *pdoflist=doflist; + return pe; } /*}}}*/ Index: /issm/trunk/src/c/objects/Loads/Numericalflux.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Numericalflux.h (revision 5910) +++ /issm/trunk/src/c/objects/Loads/Numericalflux.h (revision 5911) @@ -13,5 +13,6 @@ class Inputs; class IoModel; - +class ElementMatrix; +class ElementVector; /*}}}*/ @@ -72,10 +73,9 @@ /*}}}*/ /*Numericalflux management:{{{1*/ - void GetDofList(int** pdoflist,int approximation_enum,int setenum); void GetNormal(double* normal,double xyz_list[4][3]); - void CreateKMatrixInternal(Mat Kgg); - void CreateKMatrixBoundary(Mat Kgg); - void CreatePVectorInternal(Vec pg); - void CreatePVectorBoundary(Vec pg); + ElementMatrix* CreateKMatrixInternal(void); + ElementMatrix* CreateKMatrixBoundary(void); + ElementVector* CreatePVectorInternal(void); + ElementVector* CreatePVectorBoundary(void); /*}}}*/ Index: /issm/trunk/src/c/objects/Loads/Riftfront.cpp =================================================================== --- /issm/trunk/src/c/objects/Loads/Riftfront.cpp (revision 5910) +++ /issm/trunk/src/c/objects/Loads/Riftfront.cpp (revision 5911) @@ -378,4 +378,5 @@ void Riftfront::CreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs){ /*do nothing: */ + return; } /*}}}1*/