Changeset 25945

Timestamp:

01/19/21 03:10:32 (4 years ago)

Author:

Cheng Gong

Message:

CHG: use rotated coordinate system for Nitsches method, tested with ISMIP A-D, works for MINIcondensed and P1P1GLS. This implementation is currently only avaliable for 3D FS.

Location:

issm/trunk-jpl/src/c

Files:

• analyses/StressbalanceAnalysis.cpp (modified) (6 diffs)
• classes/Elements/Penta.cpp (modified) (3 diffs)
• cores/stressbalance_core.cpp (modified) (1 diff)

issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

@@ -5048 +5048 @@
         /*coefficient of Nitsche's method*/
         IssmDouble gamma, hx, hy, hz;
-        IssmDouble pnu, qnv;
 
         /*Get problem dimension*/
@@ -5065 +5064 @@
         IssmDouble* vdbasis = xNew<IssmDouble>(dim*vnumnodes);
         IssmDouble* pbasis  = xNew<IssmDouble>(pnumnodes);
-        IssmDouble* vdbasisn = xNew<IssmDouble>(vnumnodes);
-        IssmDouble* nsigma = xNew<IssmDouble>(numdof);
-        IssmDouble  bed_normal[3], bed_tangent[6];
 
         /*Prepare coordinate system list*/
@@ -5100 +5096 @@
         }
 
-        element->NormalBase(&bed_normal[0],xyz_list_base);
-        element->TangentBase(&bed_tangent[0], &bed_normal[0]);
         element->ElementSizes(&hx,&hy,&hz);
         element->FindParam(&gamma,FeFSNitscheGammaEnum);
@@ -5126 +5120 @@
                 element->NodalFunctionsPressure(pbasis,gauss);
                 element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
-
-                /*Frictions*/
+                
                 for(int i=0;i<vnumnodes;i++){
                         for(int j=0;j<vnumnodes;j++){
-                                for (int d=0; d<2; d++) {
-                                        for (int p=0; p<dim; p++) {
-                                                for (int q=0; q<dim; q++) {
-                                                        Ke->values[(dim*i+p)*numdof+dim*j+q] += alpha2*gauss->weight*Jdet*vbasis[i]*vbasis[j]*bed_tangent[d*3+p]*bed_tangent[d*3+q];
-                                                }
-                                        }
+                                for (int d=0; d<dim; d++) {
+                                                Ke->values[(dim*i+d)*numdof+dim*j+d] += alpha2*gauss->weight*Jdet*vbasis[i]*vbasis[j];
                                 }
                         }
                 }
                 /* -------- Nitsche terms -------- */
-                /* n\cdot\nabla\phi*/
-                for(int i=0;i<vnumnodes;i++){
-                        vdbasisn[i] = 0.0;
-                        for (int d=0; d<dim; d++) {
-                                vdbasisn[i] += bed_normal[d]*vdbasis[d*vnumnodes+i];
-                        }
-                }
                 /* Boundary terms for integration by parts.
                         The coefficient matrix of n*sigma*n-gamma*n*u is stored in the following order:
@@ -5153 +5135 @@
                 */
                 for(int i=0;i<vnumnodes;i++){
-                        for(int d=0;d<dim;d++){
-                                nsigma[d*vnumnodes+i] = 1.0 * viscosity * bed_normal[d]*vdbasisn[i];
+                        for(int j=0;j<vnumnodes;j++){
+                                /* gamma*(n*u)*(n*v) */
+                                Ke->values[(dim*i+(dim-1))*numdof+dim*j+(dim-1)] += gauss->weight*Jdet * (gamma*viscosity) * vbasis[j] * vbasis[i];
+                                /* -sigma(v)*(n*u) */
+                                Ke->values[(dim*i+(dim-1))*numdof+dim*j+(dim-1)] += (- gauss->weight*Jdet) * 2.0 * viscosity * (-vbasis[j]) * vdbasis[(dim-1)*vnumnodes+i];
+                                /* -sigma(u)*(n*v) */
+                                Ke->values[(dim*i+(dim-1))*numdof+dim*j+(dim-1)] += (- gauss->weight*Jdet) * 2.0 * viscosity * vdbasis[(dim-1)*vnumnodes+j] * (-vbasis[i]);
                         }
                 }
-                for(int i=0;i<pnumnodes;i++){
-                        nsigma[dim*vnumnodes+i] = -FSreconditioning*pbasis[i];
-                }
-
-                /*velocity components A11 */
-                for(int i=0;i<vnumnodes;i++){
-                        for (int p=0; p<dim; p++) {
-                                for(int j=0;j<vnumnodes;j++){
-                                        for (int q=0; q<dim; q++) {
-                                                /* gamma*(n*u)*(n*v) */
-                                                Ke->values[(dim*i+p)*numdof+dim*j+q] += gauss->weight*Jdet * (gamma) * bed_normal[q] * vbasis[j]*bed_normal[p] * vbasis[i];
-                                                /* sigma(u)*(n*v) */
-                                                Ke->values[(dim*i+p)*numdof+dim*j+q] += (- gauss->weight*Jdet)* nsigma[p*vnumnodes+i] * bed_normal[q] * vbasis[j];
-                                                /* sigma(v)*(n*u) */
-                                                Ke->values[(dim*i+p)*numdof+dim*j+q] += (- gauss->weight*Jdet)* bed_normal[p] * vbasis[i] * nsigma[q*vnumnodes+j];
-                                        }
-                                }
+                /* pressure x velocity  component A12, +p*(n*v) */
+                for(int k=0;k<pnumnodes;k++){
+                        for(int i=0;i<vnumnodes;i++){
+                                Ke->values[(dim*i+dim-1)*numdof+dim*vnumnodes+k] += gauss->weight*Jdet * FSreconditioning*pbasis[k] * (-vbasis[i]);
                         }
                 }
-
-                /* pressure x velocity  component A12 */
+                /* velocity x pressure component A21, +(n*u)*q */
                 for(int k=0;k<pnumnodes;k++){
-                        pnu = nsigma[dim*vnumnodes+k];
-                        for(int i=0;i<vnumnodes;i++){
-                                for (int p=0;p<dim;p++) {
-                                        Ke->values[(dim*i+p)*numdof+dim*vnumnodes+k] += gauss->weight*Jdet * pnu * (- bed_normal[p] * vbasis[i]);
-                                }
+                        for(int j=0;j<vnumnodes;j++){
+                                Ke->values[(dim*vnumnodes+k)*numdof+dim*j+dim-1] += gauss->weight*Jdet * (-vbasis[j]) * FSreconditioning*pbasis[k];
                         }
                 }
-                /* velocity x pressure component A21 */
-                for(int k=0;k<pnumnodes;k++){
-                        qnv = nsigma[dim*vnumnodes+k];
-                        for(int j=0;j<vnumnodes;j++){
-                                for (int q=0;q<dim;q++) {
-                                        Ke->values[(dim*vnumnodes+k)*numdof+dim*j+q] += gauss->weight*Jdet * ( - bed_normal[q] * vbasis[j]) * qnv;
-                                }
-                        }
-                }
-        }
-
-        /*Transform Coordinate System*/
-        element->TransformStiffnessMatrixCoord(Ke,cs_list);
+        }
 
         /*Clean up and return*/
@@ -5208 +5166 @@
         xDelete<IssmDouble>(vdbasis);
         xDelete<IssmDouble>(pbasis);
-        xDelete<IssmDouble>(vdbasisn);
-        xDelete<IssmDouble>(nsigma);
-        xDelete<int>(cs_list);
         return Ke;
 }/*}}}*/

issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

@@ -3327 +3327 @@
         int          numindices;
         int         *indices = NULL;
+        bool            isNitsche;
         IssmDouble   slopex,slopey,groundedice;
         IssmDouble   xz_plane[6];
@@ -3333 +3334 @@
         /*For FS only: we want the CS to be tangential to the bedrock*/
         this->Element::GetInputValue(&approximation,ApproximationEnum);
+        this->FindParam(&isNitsche,FlowequationIsNitscheEnum);
+
         if(!IsOnBase() || (approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum &&  approximation!=HOFSApproximationEnum)) return;
 
@@ -3357 +3360 @@
                 xz_plane[2]=slopex;   xz_plane[5]=1.;
 
-                if(groundedice>=0){
-                        if(this->nodes[indices[i]]->GetApproximation()==FSvelocityEnum){
-                                this->nodes[indices[i]]->DofInSSet(2); //vz
+                if (!isNitsche) {
+                        if(groundedice>=0){
+                                if(this->nodes[indices[i]]->GetApproximation()==FSvelocityEnum){
+                                        this->nodes[indices[i]]->DofInSSet(2); //vz
+                                }
+                                else if(this->nodes[indices[i]]->GetApproximation()==SSAFSApproximationEnum || this->nodes[indices[i]]->GetApproximation()==HOFSApproximationEnum){
+                                        this->nodes[indices[i]]->DofInSSet(4); //vz
+                                }
+                                else _error_("Flow equation approximation"<<EnumToStringx(this->nodes[indices[i]]->GetApproximation())<<" not supported yet");
                         }
-                        else if(this->nodes[indices[i]]->GetApproximation()==SSAFSApproximationEnum || this->nodes[indices[i]]->GetApproximation()==HOFSApproximationEnum){
-                                this->nodes[indices[i]]->DofInSSet(4); //vz
+                        else{
+                                if(this->nodes[indices[i]]->GetApproximation()==FSvelocityEnum){
+                                        this->nodes[indices[i]]->DofInFSet(2); //vz
+                                }
+                                else if(this->nodes[indices[i]]->GetApproximation()==SSAFSApproximationEnum || this->nodes[indices[i]]->GetApproximation()==HOFSApproximationEnum){
+                                        this->nodes[indices[i]]->DofInFSet(4); //vz
+                                }
+                                else _error_("Flow equation approximation"<<EnumToStringx(this->nodes[indices[i]]->GetApproximation())<<" not supported yet");
                         }
-                        else _error_("Flow equation approximation"<<EnumToStringx(this->nodes[indices[i]]->GetApproximation())<<" not supported yet");
-                }
-                else{
-                        if(this->nodes[indices[i]]->GetApproximation()==FSvelocityEnum){
-                                this->nodes[indices[i]]->DofInFSet(2); //vz
-                        }
-                        else if(this->nodes[indices[i]]->GetApproximation()==SSAFSApproximationEnum || this->nodes[indices[i]]->GetApproximation()==HOFSApproximationEnum){
-                                this->nodes[indices[i]]->DofInFSet(4); //vz
-                        }
-                        else _error_("Flow equation approximation"<<EnumToStringx(this->nodes[indices[i]]->GetApproximation())<<" not supported yet");
-                }
-
+                }
                 XZvectorsToCoordinateSystem(&this->nodes[indices[i]]->coord_system[0][0],&xz_plane[0]);
         }

issm/trunk-jpl/src/c/cores/stressbalance_core.cpp

@@ -47 +47 @@
                 bedslope_core(femmodel);
                 femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
-                if (!isNitsche) {
-                        ResetFSBasalBoundaryConditionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
-                }
+                ResetFSBasalBoundaryConditionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
 
                 /*We need basal melt rates for the shelf dampening*/