I managed to generate the requried mesh using gmsh as in the following code:
SetFactory("OpenCASCADE");
Rectangle (1) = {-3, -3, 0, 6, 6};
DBC() = CombinedBoundary{Surface{1};}; // Dirichlet BC
Physical Curve(1) = DBC(); // whole boundary, can use index to access different parts
Disk (2) = {0, 0, 0, 0.5};
Disk (3) = {0.8, 0, 0, 0.4};
Disk (4) = {0.3, 0.6, 0, 0.3};
v() = BooleanUnion{ Surface{2}; Delete; }{ Surface{3, 4}; Delete; }; // combine all together with conformal boundaries
w() = BooleanFragments{Surface{1}; Delete; }{Surface{v()}; Delete; };
Physical Surface(1) = w(); // whole domain, can use index to access different parts
Mesh.Algorithm = 6;
Mesh.CharacteristicLengthMin = 0.1;
Mesh.CharacteristicLengthMax = 0.1;
It has only one physical boundary, which is the boundary of the square. Then I follow this example and wrote the following code (loading .mesh or .msh file using gmshload failed, so I followed your suggestion to use DM):
load "gmsh"
load "PETSc-complex"
macro dimension()2// EOM
include "macro_ddm.idp"
// import mesh from gmsh
DM dm("3disks.msh");
if(!mpirank)
mesh Th = dm;
broadcast(processor(0), Th);
plot(Th);
Mat<complex> A;
func Pk = P2;
createMat(Th, A, Pk);
// define the bilinear form
real k = 2.0;
macro grad(u)[dx(u), dy(u)]// EOM
varf vPb(u, v) = int2d(Th)(- (grad(u)' * grad(v)) + k^2 * u * v) + on(1, u=0);
fespace Vh(Th, Pk);
A = vPb(Vh, Vh, tgv=-1);
// define the linear form
int nRhs = 1;
complex[int, int] rhs(Vh.ndof, nRhs), sol(Vh.ndof, nRhs);
for(int i = 0; i < nRhs; ++i) {
func g = (x^2+y^2<0.5^2) | ((x-0.8)^2+y^2<0.4^2) | ((x-0.3)^2+(y-0.6)^2<0.3^2); // will use multiple rhs later
varf vRhs(u, v) = int2d(Th)(g * v) + on(1, u=0);
rhs(:, i) = vRhs(0, Vh, tgv=-1);
}
// solve the system
complex[int, int] B(A.n, rhs.m), X(A.n, rhs.m);
ChangeNumbering(A, rhs, B);
set(A, sparams = "-ksp_type preonly -pc_type lu");
KSPSolve(A, B, X);
ChangeNumbering(A, sol, X, inverse = true, exchange = true);
for(int i = 0; i < rhs.m; ++i) {
Vh u;
u[] = sol(:, i).re;
plotD(Th, u, cmm="u_"+i);
}
The code runs well if i use -np 1. However, when it goes to multiple processes, the the plot only show parts of the solution. So I am afraid I can’t not output the global solution.
It looks like a mixing of some of the decomposed domains.
If I use the the built-in mesh such as
border a(t=0, 2*pi) { x=cos(t); y=sin(t); label=1; }
mesh Th;
if(!mpirank)
Th = buildmesh(a(300));
broadcast(processor(0), Th);
Then it’s ok to use multiple processes. Is there something we should take care for imported mesh when using PETSc?
