Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5737) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5738) @@ -3884,5 +3884,4 @@ } - /*Add Ke to global matrix Kgg: */ MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES); Index: /issm/trunk/src/c/objects/Elements/TriaRef.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/TriaRef.cpp (revision 5737) +++ /issm/trunk/src/c/objects/Elements/TriaRef.cpp (revision 5738) @@ -118,4 +118,44 @@ *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh3[0][i]; } +} +/*}}}*/ +/*FUNCTION TriaRef::GetSegmentBFlux{{{1*/ +void TriaRef::GetSegmentBFlux(double* B,GaussTria* gauss, int index1,int index2){ + /*Compute B matrix. B=[phi1 phi2 -phi3 -phi4] + * + * and phi1=phi3 phi2=phi4 + * + * We assume B has been allocated already, of size: 1x4 + */ + + const int numgrids=3; + double l1l3[numgrids]; + + GetNodalFunctions(&l1l3[0],gauss); + + B[0] = +l1l3[index1]; + B[1] = +l1l3[index2]; + B[2] = -l1l3[index1]; + B[3] = -l1l3[index2]; +} +/*}}}*/ +/*FUNCTION TriaRef::GetSegmentBprimeFlux{{{1*/ +void TriaRef::GetSegmentBprimeFlux(double* Bprime,GaussTria* gauss, int index1,int index2){ + /*Compute Bprime matrix. Bprime=[phi1 phi2 phi3 phi4] + * + * and phi1=phi3 phi2=phi4 + * + * We assume Bprime has been allocated already, of size: 1x4 + */ + + const int numgrids=3; + double l1l3[numgrids]; + + GetNodalFunctions(&l1l3[0],gauss); + + Bprime[0] = l1l3[index1]; + Bprime[1] = l1l3[index2]; + Bprime[2] = l1l3[index1]; + Bprime[3] = l1l3[index2]; } /*}}}*/ Index: /issm/trunk/src/c/objects/Elements/TriaRef.h =================================================================== --- /issm/trunk/src/c/objects/Elements/TriaRef.h (revision 5737) +++ /issm/trunk/src/c/objects/Elements/TriaRef.h (revision 5738) @@ -48,4 +48,6 @@ void GetNodalFunctions(double* l1l2l3,GaussTria* gauss); void GetSegmentNodalFunctions(double* l1l2l3,GaussTria* gauss, int index1,int index2); + void GetSegmentBFlux(double* B,GaussTria* gauss, int index1,int index2); + void GetSegmentBprimeFlux(double* Bprime,GaussTria* gauss, int index1,int index2); void GetNodalFunctionsDerivatives(double* l1l2l3,double* xyz_list, double* gauss); void GetNodalFunctionsDerivatives(double* l1l2l3,double* xyz_list, GaussTria* gauss); Index: /issm/trunk/src/c/objects/Loads/Numericalflux.cpp =================================================================== --- /issm/trunk/src/c/objects/Loads/Numericalflux.cpp (revision 5737) +++ /issm/trunk/src/c/objects/Loads/Numericalflux.cpp (revision 5738) @@ -3,6 +3,8 @@ */ +/*Headers:*/ +/*{{{1*/ #ifdef HAVE_CONFIG_H - #include "config.h" +#include "config.h" #else #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" @@ -15,6 +17,9 @@ #include "../../include/include.h" #include "../objects.h" - -extern int my_rank; +/*}}}*/ + +/*Load macros*/ +#define NUMVERTICES 4 +#define NDOF1 1 /*Numericalflux constructors and destructor*/ @@ -394,106 +399,59 @@ void Numericalflux::CreateKMatrixInternal(Mat Kgg){ - /* local declarations */ - int i,j,ig; - int analysis_type; - - /* node data: */ - const int numgrids=4; - const int NDOF1=1; - const int numdof=NDOF1*numgrids; - double xyz_list[numgrids][3]; - double normal[2]; - int* doflist=NULL; - - /* gaussian points: */ - int num_gauss; - double* gauss_coords =NULL; - double gauss_coord; - double* gauss_weights=NULL; - double gauss_weight; - - /* matrices: */ - double B[numgrids]; - double L[numgrids]; - double DL1,DL2; - double Ke_gg1[numdof][numdof]; - double Ke_gg2[numdof][numdof]; - double Ke_gg[numdof][numdof]={0.0}; - double Jdet; - - /*input parameters for structural analysis (diagnostic): */ - double vx,vy; - double UdotN; - double dt; - int found; - int dim; - - /*dynamic objects pointed to by hooks: */ - Input *vxaverage_input = NULL; - Input *vyaverage_input = NULL; - - /*Retrieve parameters: */ + /* constants*/ + const int numdof=NDOF1*NUMVERTICES; + + /* Intermediaries*/ + int i,j,ig,index1,index2,analysis_type; + double DL1,DL2,Jdet,dt,vx,vy,UdotN; + double xyz_list[NUMVERTICES][3]; + double normal[2]; + double B[numdof]; + double Bprime[numdof]; + double Ke_gg1[numdof][numdof]; + double Ke_gg2[numdof][numdof]; + double Ke_gg[numdof][numdof] = {0.0}; + int *doflist = NULL; + GaussTria *gauss; + + GetVerticesCoordinates(&xyz_list[0][0], nodes,NUMVERTICES); + GetDofList(&doflist); + + /*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); this->parameters->FindParam(&analysis_type,AnalysisTypeEnum); - parameters->FindParam(&dim,DimEnum); - - /*recover objects from hooks: */ - Tria* tria=(Tria*)element; - - /*recover parameters: */ - if (analysis_type==PrognosticAnalysisEnum){ - parameters->FindParam(&dt,DtEnum); - } - else if (analysis_type==BalancedthicknessAnalysisEnum || analysis_type==AdjointBalancedthicknessAnalysisEnum){ - /*No transient term is involved*/ - dt=1; - } - else{ - ISSMERROR("analysis_type %s not supported yet",EnumToString(analysis_type)); - } - - /* Get node coordinates, dof list and normal vector: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist); GetNormal(&normal[0],xyz_list); - - /*Retrieve all inputs we will be needing: */ - if(dim==2){ - vxaverage_input=tria->inputs->GetInput(VxEnum); - vyaverage_input=tria->inputs->GetInput(VyEnum); - } - - /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ - num_gauss=2; - gaussSegment(&gauss_coords, &gauss_weights,num_gauss); + switch(analysis_type){ + case PrognosticAnalysisEnum: + parameters->FindParam(&dt,DtEnum); break; + case BalancedthicknessAnalysisEnum: case AdjointBalancedthicknessAnalysisEnum: + dt=1; break;/*No transient term is involved*/ + default: + ISSMERROR("analysis_type %s not supported yet",EnumToString(analysis_type)); + } /* Start looping on the number of gaussian points: */ - for (ig=0; igGetParameterValue(&vx,gauss_coord,vxaverage_input); - this->GetParameterValue(&vy,gauss_coord,vyaverage_input); - + index1=tria->GetNodeIndex(nodes[0]); + index2=tria->GetNodeIndex(nodes[1]); + gauss=new GaussTria(index1,index2,2); + for(ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + tria->GetSegmentBFlux(&B[0],gauss,index1,index2); + tria->GetSegmentBprimeFlux(&Bprime[0],gauss,index1,index2); + + vxaverage_input->GetParameterValue(&vx,gauss); + vyaverage_input->GetParameterValue(&vy,gauss); UdotN=vx*normal[0]+vy*normal[1]; - // if (fabs(UdotN)<1.0e-9 && analysis_type==BalancedthicknessAnalysisEnum) printf("Edge number %i has a flux very small (u.n = %g ), which could lead to unaccurate results\n",id,UdotN); - - /*Get L and B: */ - GetL(&L[0],gauss_coord,numdof); - GetB(&B[0],gauss_coord); - - /*Compute DLs*/ - DL1=gauss_weight*Jdet*dt*UdotN/2; - DL2=gauss_weight*Jdet*dt*fabs(UdotN)/2; - - /* Do the triple products: */ + tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss); + DL1=gauss->weight*Jdet*dt*UdotN/2; + DL2=gauss->weight*Jdet*dt*fabs(UdotN)/2; + TripleMultiply(&B[0],1,numdof,1, &DL1,1,1,0, - &L[0],1,numdof,0, + &Bprime[0],1,numdof,0, &Ke_gg1[0][0],0); TripleMultiply(&B[0],1,numdof,1, @@ -502,17 +460,12 @@ &Ke_gg2[0][0],0); - /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ - for( i=0; i