Changeset 24714

Timestamp:

04/16/20 21:30:20 (5 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added new implementation of Dirichlet condition with zeroing rows and 1 on diagonal

File:

• issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp (modified) (3 diffs)

issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp

@@ -8 +8 @@
 #include "../modules/modules.h"
 #include "../analyses/analyses.h"
+#define USEPENALTYMETHOD false
 
 #ifdef _HAVE_PETSC_
@@ -94 +95 @@
                 MatRestoreRow(D,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
         }
-
-        /*Broadcast max(dmax)*/
-        IssmDouble dmax_all;
-        ISSM_MPI_Reduce(&dmax,&dmax_all,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&dmax_all,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-        dmax = dmax_all;
-
-        /*Penalize Dirichlet boundary*/
-        dmax = dmax * 1.e+3;
-        for(int i=0;i<femmodel->constraints->Size();i++){
-                Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
-                constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
-                if(dof!=-1){
-                        MatSetValue(LHS,dof,dof,dmax,INSERT_VALUES);
-                }
+        MatAssemblyBegin(LHS,MAT_FINAL_ASSEMBLY);
+        MatAssemblyEnd(  LHS,MAT_FINAL_ASSEMBLY);
+
+        /*Deal with Dirichlet conditions: 2 options, penalties or zeros in K matrix*/
+        if(USEPENALTYMETHOD){
+                /*Option 1: Penalty method*/
+
+                /*Broadcast max(dmax)*/
+                IssmDouble dmax_all;
+                ISSM_MPI_Reduce(&dmax,&dmax_all,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
+                ISSM_MPI_Bcast(&dmax_all,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+                dmax = dmax_all;
+
+                dmax = dmax * 1.e+3;
+                for(int i=0;i<femmodel->constraints->Size();i++){
+                        Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
+                        constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
+                        if(dof!=-1){
+                                MatSetValue(LHS,dof,dof,dmax,INSERT_VALUES);
+                        }
+                }
+        }
+        else{
+                /*Option 2: zero stiffness matrix, 1 one diagonal*/
+                int numrows = femmodel->constraints->Size();
+                int* rows = xNew<int>(numrows);
+                IssmDouble* rows_spc = xNew<IssmDouble>(numrows);
+                numrows = 0;
+
+                dmax = dmax * 1.e+3;
+                for(int i=0;i<femmodel->constraints->Size();i++){
+                        Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
+                        constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
+                        if(dof!=-1){
+                                rows[numrows] = dof;
+                                numrows++;
+                        }
+                }
+                MatZeroRows(LHS,numrows,rows,1.,NULL,NULL);
+                xDelete<int>(rows);
+                xDelete<IssmDouble>(rows_spc);
         }
         MatAssemblyBegin(LHS,MAT_FINAL_ASSEMBLY);
@@ -146 +173 @@
         VecFree(&Du);
 
-        /*Penalize Dirichlet boundary*/
-        for(int i=0;i<femmodel->constraints->Size();i++){
-                Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
-                constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
-                d = d*dmax;
-                if(dof!=-1){
-                        VecSetValues(RHS,1,&dof,(const double*)&d,INSERT_VALUES);
-                }
+        VecAssemblyBegin(RHS);
+        VecAssemblyEnd(  RHS);
+
+        /*Deal with Dirichlet conditions: 2 options, penalties or zeros in K matrix*/
+        if(USEPENALTYMETHOD){
+                /*Option 1: Penalty method*/
+                for(int i=0;i<femmodel->constraints->Size();i++){
+                        Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
+                        constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
+                        d = d*dmax;
+                        if(dof!=-1){
+                                VecSetValues(RHS,1,&dof,(const double*)&d,INSERT_VALUES);
+                        }
+                }
+        }
+        else{
+                /*Option 2: zero stiffness matrix, 1 one diagonal*/
+                int  numrows = femmodel->constraints->Size();
+                int* rows = xNew<int>(numrows);
+                IssmDouble* rows_spc = xNew<IssmDouble>(numrows);
+                numrows = 0;
+                for(int i=0;i<femmodel->constraints->Size();i++){
+                        Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
+                        constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
+                        if(dof!=-1){
+                                rows[numrows] = dof;
+                                rows_spc[numrows] = d;
+                                numrows++;
+                        }
+                }
+                VecSetValues(RHS,numrows,rows,rows_spc,INSERT_VALUES);
+                xDelete<int>(rows);
+                xDelete<IssmDouble>(rows_spc);
         }
         VecAssemblyBegin(RHS);