Skip to content

Commit

Permalink
dev merge
Browse files Browse the repository at this point in the history
  • Loading branch information
@sgarnotel
sgarnotel committed Jul 29, 2024
2 parents d6b4cd6 + 1ec138e commit 0077b0c
Show file tree
Hide file tree
Showing 5 changed files with 97 additions and 32 deletions.
6 changes: 3 additions & 3 deletions 3rdparty/getall
View file
Original file line number Diff line number Diff line change
Expand Up @@ -242,10 +242,10 @@ download(
);
download(
'PETSc',
'https://web.cels.anl.gov/projects/petsc/download/release-snapshots/petsc-3.21.2.tar.gz',
'https://web.cels.anl.gov/projects/petsc/download/release-snapshots/petsc-3.21.3.tar.gz',
'https://web.cels.anl.gov/projects/petsc/download/release-snapshots',
'petsc-3.21.2.tar.gz',
'0bb68b348861a796ddfc2e94e3fcd02d'
'petsc-3.21.3.tar.gz',
'c18cb846f97dffb9319d43e9017896b3'
);

download(
Expand Down
18 changes: 10 additions & 8 deletions examples/tutorial/tgv-test.edp
View file
Original file line number Diff line number Diff line change
@@ -1,15 +1,17 @@
load "msh3"
meshL Th = segment(3);
real [int] vtgv = [ -20,-10,-2,-1,0, tgv];
real [int] vtgv = [ -30,-20,-10,-3,-2,-1,0, tgv];
real [int] vsym = [ 0,1];
varf va(u,v)= int1d(Th)(u*v*4*9.)+on(1,2,u=0);// mul by 12 to have integer coef ..
// label : 1 on first DoF O
// 2 on last DoF 2
// tgv : -20 set to zero of row and colnum of Dof with label 1 or 2
// tgv : -10 set to zero of row of Dof with label 1 or 2
// tgv : -2 set to zero of row and colnum of Dof with label 1 or 2 and set one on diagonal term
// tgv : -1 set to zero of row of Dof with label 1 or 2 and set one on diagonal term
// tgv >=0 set to tgv value on diagonal term of Dof with label 1 or 2
// 2 on last DoF 2
// tgv : -30 set column of Dof with label 1 or 2 to zero
// tgv : -20 set row and column of Dof with label 1 or 2 to zero
// tgv : -10 set row of Dof with label 1 or 2 to zero
// tgv : -3 set column of Dof with label 1 or 2 to zero and set diagonal term to one
// tgv : -2 set row and column of Dof with label 1 or 2 to zero and set diagonal term to one
// tgv : -1 set row of Dof with label 1 or 2 to zero and set diagonal term to one
// tgv >=0 set diagonal term of Dof with label 1 or 2 to tgv value

fespace Vh(Th,P1);
int symj = 0;
Expand All @@ -23,4 +25,4 @@ for [i,tgvi:vtgv]
for [i,j,aij:A]
F(i,j) = aij;
cout << " sym= " << symj << " tgv = " << tgvi << " matrix= " << F << "\n\n\n";
}
}
83 changes: 69 additions & 14 deletions plugin/mpi/PETSc-code.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -918,12 +918,12 @@ namespace PETSc {
MatAssemblyEnd(ptA->_petsc, MAT_FINAL_ASSEMBLY);
}
if (ptParent) {
Mat** mat;
PetscInt M, N;
PetscBool assembled;
MatNestGetSubMats(ptParent->_petsc, &M, &N, &mat);
MatAssembled(ptParent->_petsc, &assembled);
if (!assembled) {
PetscInt M, N;
Mat** mat;
MatNestGetSubMats(ptParent->_petsc, &M, &N, &mat);
PetscBool assemble = PETSC_TRUE;
for (PetscInt i = 0; i < M && assemble; ++i) {
for (PetscInt j = 0; j < N && assemble; ++j) {
Expand All @@ -934,12 +934,51 @@ namespace PETSc {
}
}
}
if (assemble) {
MatAssemblyBegin(ptParent->_petsc, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(ptParent->_petsc, MAT_FINAL_ASSEMBLY);
if (assemble) assembled = PETSC_TRUE;
}
if (PetscDefined(USE_COMPLEX)) {
for (PetscInt i = 0; i < M; ++i) {
for (PetscInt j = 0; j < N; ++j) {
if (mat[i][j]) {
MatType type;
PetscBool isType[2];
MatGetType(mat[i][j], &type);
PetscStrcmp(type, MATHERMITIANTRANSPOSEVIRTUAL, isType);
if (!isType[0])
PetscStrcmp(type, MATTRANSPOSEVIRTUAL, isType + 1);
else
isType[1] = PETSC_FALSE;
if (isType[0] || isType[1]) {
Mat C, D;
if (isType[0])
MatHermitianTransposeGetMat(mat[i][j], &C);
else
MatTransposeGetMat(mat[i][j], &C);
if (C == ptA->_petsc) {
PetscReal norm;
MatDuplicate(C, MAT_COPY_VALUES, &D);
MatRealPart(D);
MatAXPY(D, -1.0, C, SAME_NONZERO_PATTERN);
MatNorm(D, NORM_INFINITY, &norm);
MatDestroy(&D);
if (norm < PETSC_MACHINE_EPSILON && isType[0]) {
MatCreateTranspose(C, &D);
MatNestSetSubMat(ptParent->_petsc, i, j, D);
MatDestroy(&D);
}
else if (norm > PETSC_MACHINE_EPSILON && isType[1]) {
MatCreateHermitianTranspose(C, &D);
MatNestSetSubMat(ptParent->_petsc, i, j, D);
MatDestroy(&D);
}
i = M, j = N;
}
}
}
}
}
}
else {
if (assembled) {
MatAssemblyBegin(ptParent->_petsc, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(ptParent->_petsc, MAT_FINAL_ASSEMBLY);
}
Expand Down Expand Up @@ -3265,7 +3304,7 @@ namespace PETSc {
}
void prepareConvert(Mat A, Mat* B) {
MatType type;
PetscBool isType;
PetscBool isType[2];
Mat** mat;
PetscInt M, N;
MatNestGetSubMats(A, &M, &N, &mat);
Expand All @@ -3275,13 +3314,23 @@ namespace PETSc {
for (PetscInt j = 0; j < N; ++j) {
if (mat[i][j]) {
MatGetType(mat[i][j], &type);
PetscStrcmp(type, MATHERMITIANTRANSPOSEVIRTUAL, &isType);
if (isType) {
PetscStrcmp(type, MATHERMITIANTRANSPOSEVIRTUAL, isType);
if (PetscDefined(USE_COMPLEX) && !isType[0])
PetscStrcmp(type, MATTRANSPOSEVIRTUAL, isType + 1);
else
isType[1] = PETSC_FALSE;
if (isType[0] || isType[1]) {
b.emplace_back(std::make_pair(std::make_pair(i, j), Mat( )));
Mat D = mat[i][j];
Mat C;
MatHermitianTransposeGetMat(D, &b.back( ).second);
MatHermitianTranspose(b.back( ).second, MAT_INITIAL_MATRIX, &C);
if (isType[0]) {
MatHermitianTransposeGetMat(D, &b.back( ).second);
MatHermitianTranspose(b.back( ).second, MAT_INITIAL_MATRIX, &C);
}
else {
MatTransposeGetMat(D, &b.back( ).second);
MatTranspose(b.back( ).second, MAT_INITIAL_MATRIX, &C);
}
PetscObjectReference((PetscObject)b.back().second);
MatNestSetSubMat(A, i, j, C);
MatDestroy(&C);
Expand Down Expand Up @@ -4763,10 +4812,16 @@ namespace PETSc {
PetscObjectTypeCompareAny((PetscObject)mat[i][j], &isType, MATMPIDENSE, MATSEQDENSE, "");
PetscInt n = 0, m = 0;
if(!isType) {
Mat C = NULL;
PetscStrcmp(type, MATHERMITIANTRANSPOSEVIRTUAL, &isType);
if(isType) {
Mat C;
if(isType)
MatHermitianTransposeGetMat(mat[i][j], &C);
else if(PetscDefined(USE_COMPLEX)) {
PetscStrcmp(type, MATTRANSPOSEVIRTUAL, &isType);
if(isType)
MatTransposeGetMat(mat[i][j], &C);
}
if(C) {
PetscObjectTypeCompareAny((PetscObject)C, &isType, MATMPIDENSE, MATSEQDENSE, "");
if(isType) MatGetSize(C, &m, &n);
type = MATHERMITIANTRANSPOSEVIRTUAL;
Expand Down
14 changes: 7 additions & 7 deletions src/femlib/HashMatrix.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1209,28 +1209,28 @@ void HashMatrix<I,R>::SetBC(char *wbc,double ttgv)
ntgv++;
if(ttgv<0)
{

if( wbc[ii] )
{
for (I k=p[ii];k<p[ii+1]; ++k)
if( j[k]==ii )
aij[k] = (std::abs(ttgv+10.0) < 1.0e-10 ? 0.0 : 1.0);
else
aij[k]=0;// put the line to Zero.
aij[k] = (ttgv < -9.0 ? 0.0 : 1.0);
else if ( std::abs(ttgv+3.0) > 1.0e-10 && std::abs(ttgv+30.0) > 1.0e-10)
aij[k]=0;// put row to zero
}

}
else
operator()(ii,ii)=ttgv;
}
if( std::abs(ttgv+2.0) < 1.0e-10 || std::abs(ttgv+20.0) < 1.0e-10) // remove also columm tgv == -2 .....
if( std::abs(ttgv+2.0) < 1.0e-10 || std::abs(ttgv+20.0) < 1.0e-10 || std::abs(ttgv+3.0) < 1.0e-10 || std::abs(ttgv+30.0) < 1.0e-10) // remove also columm tgv == -2,-3 .....
{
CSC();
for(I jj=0; jj< this->n; ++jj)
if( wbc[jj] ) {
for (I k=p[jj];k<p[jj+1]; ++k)
if( i[k]!=jj || std::abs(ttgv+20.0) < 1.0e-10)
aij[k]=0;//
if( i[k]!=jj || ttgv < -19.0)
aij[k]=0;// put column to zero
}
}

Expand Down
8 changes: 8 additions & 0 deletions src/fflib/problem.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9904,6 +9904,8 @@ bool AssembleVarForm(Stack stack,const MMesh & Th,const FESpace1 & Uh,const FESp
tgv = -10;
else if (std::abs(tgv+2.0) < 1.0e-10)
tgv = -20;
else if (std::abs(tgv+3.0) < 1.0e-10)
tgv = -30;
}

KN<R> gg(buf);
Expand Down Expand Up @@ -10061,6 +10063,8 @@ bool AssembleVarForm(Stack stack,const MMesh & Th,const FESpace1 & Uh,const FESp
tgv = -10;
else if (std::abs(tgv+2.0) < 1.0e-10)
tgv = -20;
else if (std::abs(tgv+3.0) < 1.0e-10)
tgv = -30;
}

int Nbcomp=Vh.N;
Expand Down Expand Up @@ -10218,6 +10222,8 @@ bool AssembleVarForm(Stack stack,const MMesh & Th,const FESpace1 & Uh,const FESp
tgv = -10;
else if (std::abs(tgv+2.0) < 1.0e-10)
tgv = -20;
else if (std::abs(tgv+3.0) < 1.0e-10)
tgv = -30;
}

int Nbcomp=Vh.N;
Expand Down Expand Up @@ -10383,6 +10389,8 @@ bool AssembleVarForm(Stack stack,const MMesh & Th,const FESpace1 & Uh,const FESp
tgv = -10;
else if (std::abs(tgv+2.0) < 1.0e-10)
tgv = -20;
else if (std::abs(tgv+3.0) < 1.0e-10)
tgv = -30;
}

int Nbcomp=Vh.N;
Expand Down

0 comments on commit 0077b0c

Please sign in to comment.