Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 24130) +++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 24131) @@ -225,75 +225,4 @@ } -} -/*}}}*/ -void TriaRef::GetSegmentBFlux(IssmDouble* B,Gauss* gauss, int index1,int index2,int finiteelement){/*{{{*/ - /*Compute B matrix. B=[phi1 phi2 -phi3 -phi4] - * - * and phi1=phi3 phi2=phi4 - * - * We assume B has been allocated already, of size: 1x4 - */ - - /*Fetch number of nodes for this finite element*/ - int numnodes = this->NumberofNodes(finiteelement); - - /*Get nodal functions*/ - IssmDouble* basis=xNew(numnodes); - GetNodalFunctions(basis,gauss,finiteelement); - - /*Build B for this segment*/ - switch(finiteelement){ - case P0DGEnum: - B[0] = +basis[0]; - B[1] = -basis[0]; - break; - case P1DGEnum: - B[0] = +basis[index1]; - B[1] = +basis[index2]; - B[2] = -basis[index1]; - B[3] = -basis[index2]; - break; - default: - _error_("not supported yet"); - } - - /*Clean-up*/ - xDelete(basis); -} -/*}}}*/ -void TriaRef::GetSegmentBprimeFlux(IssmDouble* Bprime,Gauss* gauss, int index1,int index2,int finiteelement){/*{{{*/ - /*Compute Bprime matrix. Bprime=[phi1 phi2 phi3 phi4] - * - * and phi1=phi3 phi2=phi4 - * - * We assume Bprime has been allocated already, of size: 1x4 - */ - - /*Fetch number of nodes for this finite element*/ - int numnodes = this->NumberofNodes(finiteelement); - - /*Get nodal functions*/ - IssmDouble* basis=xNew(numnodes); - GetNodalFunctions(basis,gauss,finiteelement); - - /*Build B'*/ - /*Build B for this segment*/ - switch(finiteelement){ - case P0DGEnum: - Bprime[0] = basis[0]; - Bprime[1] = basis[0]; - break; - case P1DGEnum: - Bprime[0] = basis[index1]; - Bprime[1] = basis[index2]; - Bprime[2] = basis[index1]; - Bprime[3] = basis[index2]; - break; - default: - _error_("not supported yet"); - } - - /*Clean-up*/ - xDelete(basis); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h (revision 24130) +++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h (revision 24131) @@ -25,6 +25,4 @@ void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss,int finiteelement); void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss,int finiteelement); - void GetSegmentBFlux(IssmDouble* B,Gauss* gauss, int index1,int index2,int finiteelement); - void GetSegmentBprimeFlux(IssmDouble* Bprime,Gauss* gauss, int index1,int index2,int finiteelement); void GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss); void GetSegmentNodalFunctions(IssmDouble* basis,Gauss* gauss, int index1,int index2,int finiteelement); Index: /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp (revision 24130) +++ /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp (revision 24131) @@ -574,15 +574,18 @@ /* Intermediaries*/ - IssmDouble DL1,DL2,Jdet,vx,vy,UdotN; + IssmDouble A1,A2,Jdet,vx,vy,UdotN; IssmDouble xyz_list[NUMVERTICES][3]; IssmDouble normal[2]; /*Fetch number of nodes for this flux*/ - int numnodes = this->GetNumberOfNodes(); + int numnodes = this->GetNumberOfNodes(); + int numnodes_plus = this->GetNumberOfNodesOneSide(); + int numnodes_minus = numnodes_plus; /*For now we are not doing p-adaptive DG*/ + _assert_(numnodes==numnodes_plus+numnodes_minus); /*Initialize variables*/ - ElementMatrix *Ke = new ElementMatrix(nodes,numnodes,this->parameters); - IssmDouble *B = xNew(numnodes); - IssmDouble *Bprime = xNew(numnodes); + ElementMatrix *Ke = new ElementMatrix(nodes,numnodes,this->parameters); + IssmDouble *basis_plus = xNew(numnodes_plus); + IssmDouble *basis_minus = xNew(numnodes_minus); /*Retrieve all inputs and parameters*/ @@ -600,6 +603,6 @@ gauss->GaussPoint(ig); - tria->GetSegmentBFlux(&B[0],gauss,index1,index2,tria->FiniteElement()); - tria->GetSegmentBprimeFlux(&Bprime[0],gauss,index1,index2,tria->FiniteElement()); + tria->GetSegmentNodalFunctions(&basis_plus[0] ,gauss,index1,index2,tria->FiniteElement()); + tria->GetSegmentNodalFunctions(&basis_minus[0],gauss,index1,index2,tria->FiniteElement()); vxaverage_input->GetInputValue(&vx,gauss); @@ -607,145 +610,6 @@ UdotN=vx*normal[0]+vy*normal[1]; tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss); - DL1=gauss->weight*Jdet*UdotN/2; - DL2=gauss->weight*Jdet*fabs(UdotN)/2; - - for(int i=0;ivalues[i*numnodes+j]+=DL1*B[i]*Bprime[j]; - Ke->values[i*numnodes+j]+=DL2*B[i]*B[j]; - } - } - } - - /*Clean up and return*/ - xDelete(B); - xDelete(Bprime); - delete gauss; - return Ke; -} -/*}}}*/ -ElementMatrix* Numericalflux::CreateKMatrixMasstransport(void){/*{{{*/ - - switch(this->flux_type){ - case InternalEnum: - return CreateKMatrixMasstransportInternal(); - case BoundaryEnum: - return CreateKMatrixMasstransportBoundary(); - default: - _error_("type not supported yet"); - } -} -/*}}}*/ -ElementMatrix* Numericalflux::CreateKMatrixMasstransportBoundary(void){/*{{{*/ - - /*Initialize Element matrix and return if necessary*/ - Tria* tria=(Tria*)element; - if(!tria->IsIceInElement()) return NULL; - - /* Intermediaries*/ - IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble normal[2]; - - /*Retrieve all inputs and parameters*/ - GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - IssmDouble dt = parameters->FindParam(TimesteppingTimeStepEnum); - Input* vxaverage_input=tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input); - Input* vyaverage_input=tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input); - GetNormal(&normal[0],xyz_list); - - /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/ - int index1=tria->GetVertexIndex(vertices[0]); - int index2=tria->GetVertexIndex(vertices[1]); - - GaussTria* gauss=new GaussTria(); - gauss->GaussEdgeCenter(index1,index2); - vxaverage_input->GetInputValue(&mean_vx,gauss); - vyaverage_input->GetInputValue(&mean_vy,gauss); - delete gauss; - - IssmDouble UdotN=mean_vx*normal[0]+mean_vy*normal[1]; - if(UdotN<=0){ - return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/ - } - - /*Initialize Element vector and other vectors*/ - int numnodes = this->GetNumberOfNodes(); - ElementMatrix *Ke = new ElementMatrix(nodes,numnodes,this->parameters); - IssmDouble *basis = xNew(numnodes); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(index1,index2,2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement()); - - vxaverage_input->GetInputValue(&vx,gauss); - vyaverage_input->GetInputValue(&vy,gauss); - UdotN=vx*normal[0]+vy*normal[1]; - tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss); - DL=gauss->weight*Jdet*dt*UdotN; - - for(int i=0;ivalues[i*numnodes+j]+=DL*basis[i]*basis[j]; - } - } - } - - /*Clean up and return*/ - xDelete(basis); - delete gauss; - return Ke; -} -/*}}}*/ -ElementMatrix* Numericalflux::CreateKMatrixMasstransportInternal(void){/*{{{*/ - - /*Initialize Element matrix and return if necessary*/ - Tria* tria=(Tria*)element; - if(!tria->IsIceInElement()) return NULL; - - /* Intermediaries*/ - IssmDouble A1,A2,Jdet,vx,vy,UdotN; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble normal[2]; - - /*Fetch number of nodes for this flux*/ - int numnodes = this->GetNumberOfNodes(); - int numnodes_plus = this->GetNumberOfNodesOneSide(); - int numnodes_minus = numnodes_plus; /*For now we are not doing p-adaptive DG*/ - _assert_(numnodes==numnodes_plus+numnodes_minus); - - /*Initialize variables*/ - ElementMatrix *Ke = new ElementMatrix(nodes,numnodes,this->parameters); - IssmDouble *basis_plus = xNew(numnodes_plus); - IssmDouble *basis_minus = xNew(numnodes_minus); - - /*Retrieve all inputs and parameters*/ - GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - IssmDouble dt = parameters->FindParam(TimesteppingTimeStepEnum); - Input* vxaverage_input=tria->inputs->GetInput(VxEnum); - Input* vyaverage_input=tria->inputs->GetInput(VyEnum); - GetNormal(&normal[0],xyz_list); - - /* Start looping on the number of gaussian points: */ - int index1=tria->GetVertexIndex(vertices[0]); - int index2=tria->GetVertexIndex(vertices[1]); - GaussTria* gauss=new GaussTria(index1,index2,2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - tria->GetSegmentNodalFunctions(&basis_plus[0] ,gauss,index1,index2,tria->FiniteElement()); - tria->GetSegmentNodalFunctions(&basis_minus[0],gauss,index1,index2,tria->FiniteElement()); - - vxaverage_input->GetInputValue(&vx,gauss); - vyaverage_input->GetInputValue(&vy,gauss); - UdotN=vx*normal[0]+vy*normal[1]; - tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss); - A1=gauss->weight*Jdet*dt*UdotN/2; - A2=gauss->weight*Jdet*dt*fabs(UdotN)/2; + A1=gauss->weight*Jdet*UdotN/2; + A2=gauss->weight*Jdet*fabs(UdotN)/2; /*Term 1 (numerical flux): {Hv}.[[phi]] = 0.5(H+v+ + H-v-)(phi+n+ + phi-n-) @@ -798,4 +662,177 @@ } /*}}}*/ +ElementMatrix* Numericalflux::CreateKMatrixMasstransport(void){/*{{{*/ + + switch(this->flux_type){ + case InternalEnum: + return CreateKMatrixMasstransportInternal(); + case BoundaryEnum: + return CreateKMatrixMasstransportBoundary(); + default: + _error_("type not supported yet"); + } +} +/*}}}*/ +ElementMatrix* Numericalflux::CreateKMatrixMasstransportBoundary(void){/*{{{*/ + + /*Initialize Element matrix and return if necessary*/ + Tria* tria=(Tria*)element; + if(!tria->IsIceInElement()) return NULL; + + /* Intermediaries*/ + IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble normal[2]; + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + IssmDouble dt = parameters->FindParam(TimesteppingTimeStepEnum); + Input* vxaverage_input=tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input); + Input* vyaverage_input=tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input); + GetNormal(&normal[0],xyz_list); + + /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/ + int index1=tria->GetVertexIndex(vertices[0]); + int index2=tria->GetVertexIndex(vertices[1]); + + GaussTria* gauss=new GaussTria(); + gauss->GaussEdgeCenter(index1,index2); + vxaverage_input->GetInputValue(&mean_vx,gauss); + vyaverage_input->GetInputValue(&mean_vy,gauss); + delete gauss; + + IssmDouble UdotN=mean_vx*normal[0]+mean_vy*normal[1]; + if(UdotN<=0){ + return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/ + } + + /*Initialize Element vector and other vectors*/ + int numnodes = this->GetNumberOfNodes(); + ElementMatrix *Ke = new ElementMatrix(nodes,numnodes,this->parameters); + IssmDouble *basis = xNew(numnodes); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(index1,index2,2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement()); + + vxaverage_input->GetInputValue(&vx,gauss); + vyaverage_input->GetInputValue(&vy,gauss); + UdotN=vx*normal[0]+vy*normal[1]; + tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss); + DL=gauss->weight*Jdet*dt*UdotN; + + for(int i=0;ivalues[i*numnodes+j]+=DL*basis[i]*basis[j]; + } + } + } + + /*Clean up and return*/ + xDelete(basis); + delete gauss; + return Ke; +} +/*}}}*/ +ElementMatrix* Numericalflux::CreateKMatrixMasstransportInternal(void){/*{{{*/ + + /*Initialize Element matrix and return if necessary*/ + Tria* tria=(Tria*)element; + if(!tria->IsIceInElement()) return NULL; + + /* Intermediaries*/ + IssmDouble A1,A2,Jdet,vx,vy,UdotN; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble normal[2]; + + /*Fetch number of nodes for this flux*/ + int numnodes = this->GetNumberOfNodes(); + int numnodes_plus = this->GetNumberOfNodesOneSide(); + int numnodes_minus = numnodes_plus; /*For now we are not doing p-adaptive DG*/ + _assert_(numnodes==numnodes_plus+numnodes_minus); + + /*Initialize variables*/ + ElementMatrix *Ke = new ElementMatrix(nodes,numnodes,this->parameters); + IssmDouble *basis_plus = xNew(numnodes_plus); + IssmDouble *basis_minus = xNew(numnodes_minus); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + IssmDouble dt = parameters->FindParam(TimesteppingTimeStepEnum); + Input* vxaverage_input=tria->inputs->GetInput(VxEnum); + Input* vyaverage_input=tria->inputs->GetInput(VyEnum); + GetNormal(&normal[0],xyz_list); + + /* Start looping on the number of gaussian points: */ + int index1=tria->GetVertexIndex(vertices[0]); + int index2=tria->GetVertexIndex(vertices[1]); + GaussTria* gauss=new GaussTria(index1,index2,2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + tria->GetSegmentNodalFunctions(&basis_plus[0] ,gauss,index1,index2,tria->FiniteElement()); + tria->GetSegmentNodalFunctions(&basis_minus[0],gauss,index1,index2,tria->FiniteElement()); + + vxaverage_input->GetInputValue(&vx,gauss); + vyaverage_input->GetInputValue(&vy,gauss); + UdotN=vx*normal[0]+vy*normal[1]; + tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss); + A1=gauss->weight*Jdet*dt*UdotN/2; + A2=gauss->weight*Jdet*dt*fabs(UdotN)/2; + + /*Term 1 (numerical flux): {Hv}.[[phi]] = 0.5(H+v+ + H-v-)(phi+n+ + phi-n-) + * = v.n/2 (H+phi+ + H-phi+ -H+phi- -H-phi-) + * = v.n/2 (H+phi+ + H-phi+ -H+phi- -H-phi-) + * + *Term 2 (stabilization) |v.n|/2 [[H]].[[phi]] = |v.n|/2 (H+n+ + H-n-)(phi+n+ + phi-n-) + * = |v.n|/2 (H+phi+ -H-phi+ -H+phi- +H-phi-) + * | A++ | A+- | + * K = |-----------| + * | A-+ | A-- | + * + *These 4 terms for each expressions are added independently*/ + + /*First term A++*/ + for(int i=0;ivalues[i*numnodes+j] += A1*(basis_plus[j]*basis_plus[i]); + Ke->values[i*numnodes+j] += A2*(basis_plus[j]*basis_plus[i]); + } + } + /*Second term A+-*/ + for(int i=0;ivalues[i*numnodes+numnodes_plus+j] += A1*(basis_minus[j]*basis_plus[i]); + Ke->values[i*numnodes+numnodes_plus+j] += -A2*(basis_minus[j]*basis_plus[i]); + } + } + /*Third term A-+*/ + for(int i=0;ivalues[(numnodes_plus+i)*numnodes+j] += -A1*(basis_plus[j]*basis_minus[i]); + Ke->values[(numnodes_plus+i)*numnodes+j] += -A2*(basis_plus[j]*basis_minus[i]); + } + } + /*Fourth term A-+*/ + for(int i=0;ivalues[(numnodes_plus+i)*numnodes+numnodes_plus+j] += -A1*(basis_minus[j]*basis_minus[i]); + Ke->values[(numnodes_plus+i)*numnodes+numnodes_plus+j] += A2*(basis_minus[j]*basis_minus[i]); + } + } + } + + /*Clean up and return*/ + xDelete(basis_plus); + xDelete(basis_minus); + delete gauss; + return Ke; +} +/*}}}*/ ElementVector* Numericalflux::CreatePVectorAdjointBalancethickness(void){/*{{{*/