Changeset 18349

Timestamp:

08/08/14 13:33:36 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: added CreateLHS

File:

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

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

@@ -161 +161 @@
 
         /*Create LHS: [ML − theta*detlat *(K+D)^n+1]*/
-        MatDuplicate(K_petsc,MAT_SHARE_NONZERO_PATTERN,&LHS);
-        MatGetOwnershipRange(K_petsc,&rstart,&rend);
-        for(int row=rstart; row<rend; row++){
-                MatGetRow(K_petsc,row,&ncols, (const int**)&cols, (const double**)&vals);
-                MatGetRow(D_petsc,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
-                _assert_(ncols==ncols2);
-                for(int j=0; j<ncols; j++) {
-                        _assert_(cols[j]==cols2[j]);
-                        d = -theta*deltat*(vals[j] + vals2[j]);
-                        if(cols[j]==row){
-                                VecGetValues(Ml_petsc,1,(const int*)&cols[j],&mi);
-                                d += mi;
-                        }
-                        if(fabs(d)>dmax) dmax = fabs(d);
-                        MatSetValues(LHS,1,&row,1,&cols[j],(const double*)&d,INSERT_VALUES);
-                }
-                MatRestoreRow(K_petsc,row,&ncols, (const int**)&cols, (const double**)&vals);
-                MatRestoreRow(D_petsc,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
-        }
-
-        /*Penalize Dirichlet boundary*/
-        dmax = dmax * 1.e+3;
-        for(int i=0;i<femmodel->constraints->Size();i++){
-                Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
-                if(constraint->InAnalysis(analysis_type)){
-                        constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
-                        if(dof!=-1){
-                                MatSetValues(LHS,1,&dof,1,&dof,(const double*)&dmax,INSERT_VALUES);
-                        }
-                }
-        }
-        MatAssemblyBegin(LHS,MAT_FINAL_ASSEMBLY);
-        MatAssemblyEnd(  LHS,MAT_FINAL_ASSEMBLY);
+        CreateLHS(&LHS,&dmax,K_petsc,D_petsc,Ml_petsc,theta,deltat,femmodel,configuration_type);
 
         /*Create RHS: [ML + (1 − theta) deltaT L^n] u^n */
@@ -265 +233 @@
         VecDuplicate(u,&Ri_plus);
         VecDuplicate(u,&Ri_minus);
+        MatGetOwnershipRange(K_petsc,&rstart,&rend);
         for(int row=rstart; row<rend; row++){
                 double Pi_plus  = 0.;