Hi Ohi,
I haven't had a chance to look at what you sent in detail but I suggest
that you
could start by using Overture/tests/tcm3.C to solve Laplace equation with
Dirichlet or Neumann BC's on the grid you created. You could then make
changes to get the BC's that you want.
Regards
Bill
On Wed, Jun 14, 2017 at 6:14 PM, Ohiremen L Dibua <odibua@xxxxxxxxxxxx>
wrote:
Good evening Bill,
Thank you for being so good about providing this open-source software, and
this forum to ask questions about it. I recently re-meshed a prior set-up
by creating more grids, and implementing what looks to me like proper code
for solving it. I am solving a Laplacian equation. My set up has 9 grids. I
get an error saying
Oges:generateMatrix:ERROR matrix has a row with all zero coefficients
The offending equation is (i1,i2,i3,n,grid)=(10,0,0,0,0).
Below is my code up until the point where an error is thrown, and below
that is the cmd I use to generate the relevant composite grid.
*LAPLACIAN CODE*
// create and read in a CompositeGrid
CompositeGrid cg;
getFromADataBase(cg,nameOfOGFile);
cg.update();
cg.update(MappedGrid::THEvertex | MappedGrid::THEcenter | MappedGrid::
THEvertexBoundaryNormal);
Ogshow show(nameOfShowFile);
show.saveGeneralComment("Solution to air comb-drive");
show.saveGeneralComment("Test file written on February");
// make a grid function to hold coefficients
Range all;
int stencilSize=int( pow(3,cg.numberOfDimensions())+1.5); //Add 1 for
interpolation equations
realCompositeGridFunction coeff(cg,stencilSize,all,all,all);
realCompositeGridFunction q(cg,all,all,all,2);
coeff.setIsACoefficientMatrix(TRUE,stencilSize);
q=0.;
// create grid functions:
realCompositeGridFunction u, f;
u.link(q,Range(0,0));
f.link(q,Range(1,1));
q.setName("q");
q.setName("u",0);
//q.setName("f",1);
cout << "Number of Grids " << cg.numberOfComponentGrids() << endl;
char buffer[80];
Index I1,I2,I3,I1F,I2F,I3F;
Index Ib1,Ib2,Ib3,Ib1F,Ib2F,Ib3F;
Index Ig1,Ig2,Ig3;
Index Ip1,Ip2,Ip3;
CompositeGridOperators op(cg); // create some differential operators
op.setStencilSize(stencilSize);
coeff.setOperators(op);
cout << "Here " << endl;
coeff=op.laplacianCoefficients(); // get the coefficients for the Laplace
operator
// make shorter names for readability
BCTypes::BCNames dirichlet = BCTypes::dirichlet,
neumann = BCTypes::neumann,
extrapolate = BCTypes::extrapolate,
evenSymmetry = BCTypes::evenSymmetry,
allBoundaries = BCTypes::allBoundaries;
cout << "Here 2 " << endl;
const int ltDomain=1, rtDomain=2, lowDomain=3, upDomain=4;
int grid;
grid=0;
MappedGrid & mg0 = cg[grid];
mg0.update(MappedGrid::THEvertex | MappedGrid::THEcenter | MappedGrid::
THEvertexBoundaryNormal);
mg0.setBoundaryCondition(Start,axis1,ltDomain);
mg0.setBoundaryCondition(End,axis1,rtDomain);
mg0.setBoundaryCondition(Start,axis2,lowDomain);
mg0.setBoundaryCondition(End,axis2,upDomain);
coeff[grid].applyBoundaryConditionCoefficients(0,0,dirichlet,ltDomain);
coeff[grid].applyBoundaryConditionCoefficients(0,0,extrapolate,ltDomain);
coeff[grid].applyBoundaryConditionCoefficients(0,0,dirichlet,lowDomain);
coeff[grid].applyBoundaryConditionCoefficients(0,0,extrapolate,lowDomain);
coeff[grid].applyBoundaryConditionCoefficients(0,0,neumann,upDomain);
coeff[grid].finishBoundaryConditions();
getIndex(mg0.indexRange(),I1,I2,I3);
f[grid](I1,I2,I3)=0;
*CMD FOR COMPOSITE GRID*
create mappings
#Create Domain 1 ###############
rectangle
mappingName
domain1
set corners
-12.5 2.5 0. 1.5
lines
10 10
boundary conditions
1 0 1 1
exit
#Create Domain 4 ###############
rectangle
mappingName
domain4
set corners
2.5 7.5 0 0.5
#-12.5 7.5 0 0.5
lines
10 10
boundary conditions
0 0 1 0
exit
#Create Domain 5 ###############
rectangle
mappingName
domain5
set corners
7.5 12.5 0 1
lines
10 10
boundary conditions
0 0 1 1
exit
#Create Domain 3 ###############
spline
enter spline points
5
7. 1.5
7.125 1.169
7.25 1.067
7.375 1.016
7.5 1.
lines
10
mappingName
topSpline
exit
spline
enter spline points
7
7.5 1.
7.5 0.9
7.5 0.85
7.5 0.8
7.5 0.7
7.5 0.6
7.5 0.5
lines
10
mappingName
rightSpline
exit
spline
enter spline points
9
6.5 1.5
6.55 1.188
6.6 1.064
6.7 0.9
6.8 0.789
7. 0.634
7.2 0.546
7.4 0.505
7.5 0.5
lines
10
mappingName
bottomSpline
exit
spline
enter spline points
6
6.5 1.5
6.55 1.5
6.6 1.5
6.7 1.5
6.9 1.5
7. 1.5
lines
10
mappingName
leftSpline
exit
tfi
mappingName
domain3
choose bottom
bottomSpline
choose top
topSpline
choose left
leftSpline
choose right
rightSpline
lines
10 10
boundary conditions
1 0 0 1
exit
check
domain3
#Create Domain 2 ###############
spline
enter spline points
9
6.5 1.5
6.55 1.188
6.6 1.064
6.7 0.9
6.8 0.789
7. 0.634
7.2 0.546
7.4 0.505
7.5 0.5
lines
10
mappingName
rightSplineD2
exit
spline
enter spline points
8
2.5 1.5
2.8 1.5
3.1 1.5
3.4 1.5
3.8 1.5
4.2 1.5
5.6 1.5
6.5 1.5
lines
10
mappingName
topSplineD2
exit
spline
enter spline points
9
2.5 0.5
2.8 0.5
3.1 0.5
3.4 0.5
4 0.5
5.5 0.5
6 0.5
6.6 0.5
7.5 0.5
lines
10
mappingName
bottomSplineD2
exit
spline
enter spline points
9
2.5 1.5
2.5 1.3
2.5 1.2
2.5 1.0
2.5 0.8
2.5 0.7
2.5 0.6
2.5 0.55
2.5 0.5
lines
10
mappingName
leftSplineD2
exit
tfi
mappingName
domain2
choose bottom
bottomSplineD2
choose top
topSplineD2
choose left
leftSplineD2
choose right
rightSplineD2
lines
10 10
boundary conditions
0 0 1 0
exit
check
domain2
#Create Domain 6 ###############
spline
enter spline points
5
13. -0.5
12.875 -0.169
12.75 -0.067
12.625 -0.016
12.5 0.
lines
10
mappingName
bottomSplineRight
exit
spline
enter spline points
7
12.5 0.5
12.5 0.45
12.5 0.4
12.5 0.3
12.5 0.2
12.5 0.1
12.5 0.
lines
10
mappingName
leftSplineRight
exit
spline
enter spline points
9
13.5 -0.5
13.45 -0.188
13.4 -0.0641
13.3 0.1
13.2 0.214
13 0.366
12.8 0.454
12.6 0.495
12.5 0.5
lines
10
mappingName
topSplineRight
exit
spline
enter spline points
6
15.5 -0.5
15.55 -0.5
15.6 -0.5
15.8 -0.5
15.9 -0.5
16 -0.5
lines
10
mappingName
rightSplineRight
exit
tfi
mappingName
domain6
choose bottom
bottomSplineRight
choose top
topSplineRight
choose left
leftSplineRight
choose right
rightSplineRight
lines
10 10
boundary conditions
0 1 0 1
exit
check
domain6
#Create Domain 7 ###############
spline
enter spline points
9
13.5 -0.5
13.45 -0.188
13.4 -0.0641
13.3 0.1
13.2 0.214
13 0.366
12.8 0.454
12.6 0.495
12.5 0.5
lines
10
mappingName
leftSplineRightD2
exit
spline
enter spline points
8
13.5 -0.5
14 -0.5
14.5 -0.5
15 -0.5
15.5 -0.5
16 -0.5
16.5 -0.5
17.5 -0.5
line
10
mappingName
bottomSplineRightD2
exit
spline
enter spline points
9
12.5 0.5
13.5 0.5
14 0.5
14.5 0.5
15 0.5
15.5 0.5
16 0.5
16.5 0.5
17.5 0.5
lines
10
mappingName
topSplineRightD2
exit
spline
enter spline points
9
17.5 0.5
17.5 0.3
17.5 0.2
17.5 0.
17.5 -0.1
17.5 -0.3
17.5 -0.4
17.5 -0.45
17.5 -0.5
lines
10
mappingName
rightSplineRightD2
exit
tfi
mappingName
domain7
choose bottom
bottomSplineRightD2
choose top
topSplineRightD2
choose left
leftSplineRightD2
choose right
rightSplineRightD2
lines
10 10
boundary conditions
0 0 1 0
exit
check
domain7
#Create Domain 8 ###############
rectangle
mappingName
domain8
set corners
17.5 32.5 -0.5 1.
#17.5 32.5 -0.5 0.5
lines
10 10
boundary conditions
0 1 1 1
exit
#Create Domain 9 ###############
rectangle
mappingName
domain9
set corners
12.5 17.5 0.5 1.
#12.5 32.5 0.5 1.
lines
10 10
boundary conditions
0 0 0 1
exit
exit
generate an overlapping grid
domain1
domain2
domain3
domain4
domain5
domain6
domain7
domain8
domain9
done
change parameters
prevent hole cutting
#domain1
#all
#domain2
#all
#domain3
#all
#domain4
#all
domain5
all
#domain6
#all
#domain7
#all
#domain8
#all
#domain9
#all
done
ghost points
all
1 1 1 1
exit
compute overlap
exit
save an overlapping grid
combGridNewMeshOverlap5.hdf
combGridNewMeshOverlap5
exit
Thank you for your time,
Ohi
