Changeset 24130

Timestamp:

08/30/19 11:51:24 (6 years ago)

Author:

Mathieu Morlighem

Message:

CHG: better way of computing stiffness with DG

Location:

issm/trunk-jpl/src/c/classes/Loads

Files:

• Numericalflux.cpp (modified) (4 diffs)
• Numericalflux.h (modified) (1 diff)

issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp

@@ -351 +351 @@
                                 _error_("Numericalflux type " << EnumToStringx(this->flux_type) << " not supported yet");
                 }
+        }
+        else{
+                _error_("Numericalflux " << EnumToStringx(this->flux_degree) << " not supported yet");
+        }
+
+}
+/*}}}*/
+int   Numericalflux::GetNumberOfNodesOneSide(void){/*{{{*/
+
+        if(this->flux_degree==P0DGEnum){
+                return 1;
+        }
+        else if(this->flux_degree==P1DGEnum){
+                return 2;
         }
         else{
@@ -695 +709 @@
 
         /* 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<IssmDouble>(numnodes);
-        IssmDouble    *Bprime = xNew<IssmDouble>(numnodes);
+        ElementMatrix *Ke = new ElementMatrix(nodes,numnodes,this->parameters);
+        IssmDouble    *basis_plus  = xNew<IssmDouble>(numnodes_plus);
+        IssmDouble    *basis_minus = xNew<IssmDouble>(numnodes_minus);
 
         /*Retrieve all inputs and parameters*/
@@ -722 +739 @@
                 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);
@@ -729 +746 @@
                 UdotN=vx*normal[0]+vy*normal[1];
                 tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
-                DL1=gauss->weight*Jdet*dt*UdotN/2;
-                DL2=gauss->weight*Jdet*dt*fabs(UdotN)/2;
-
-                for(int i=0;i<numnodes;i++){
-                        for(int j=0;j<numnodes;j++){
-                                Ke->values[i*numnodes+j]+=DL1*B[i]*Bprime[j];
-                                Ke->values[i*numnodes+j]+=DL2*B[i]*B[j];
+                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;i<numnodes_plus;i++){
+                        for(int j=0;j<numnodes_plus;j++){
+                                Ke->values[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;i<numnodes_plus;i++){
+                        for(int j=0;j<numnodes_minus;j++){
+                                Ke->values[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;i<numnodes_minus;i++){
+                        for(int j=0;j<numnodes_plus;j++){
+                                Ke->values[(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;i<numnodes_minus;i++){
+                        for(int j=0;j<numnodes_minus;j++){
+                                Ke->values[(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<IssmDouble>(B);
-   xDelete<IssmDouble>(Bprime);
+   xDelete<IssmDouble>(basis_plus);
+   xDelete<IssmDouble>(basis_minus);
         delete gauss;
         return Ke;

issm/trunk-jpl/src/c/classes/Loads/Numericalflux.h

@@ -64 +64 @@
                 void GetNodesSidList(int* sidlist);
                 int  GetNumberOfNodes(void);
+                int  GetNumberOfNodesOneSide(void);
                 bool IsPenalty(void);
                 void PenaltyCreateJacobianMatrix(Matrix<IssmDouble>* Jff,IssmDouble kmax){_error_("Not implemented yet");};