Am 29.03.2012 16:27, schrieb Daniel Wheeler: > On Thu, Mar 29, 2012 at 3:06 AM, Andreas Hasenkopf > <[email protected] <mailto:[email protected]>> wrote: > > Dear FiPy-users, > in my latest simulation on a mesh containing multiple volumes created by > extruding a surface in gmsh I noticed a rather strong discontinuity of > the solution (in this case solution to laplace-equation in 3D). > As you can see in the attached image (scalar cut plane, z plane is in > the image plane) it appears that the PDE is solved for each volume in > the grid independently with no regard for the solution in adjacent > volumes. > I'm wondering now whether this is a problem due to my "poor resolution" > or whether I have to generate a mesh containing only one volume element? > > > It looks as though the meshes are disconnected or the faces don't align. > Can you check the connectivity in FiPy? Maybe I can take a look if you > send the script that you are running and the geo file that generated the > mesh. > > Cheers. > > -- > Daniel Wheeler > > > _______________________________________________ > fipy mailing list > [email protected] > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Hi Daniel, thanks for having a look at the problematic geo file, too. In the meanwhile I was redrawing the geometry in another file and noticed that one ruled surface had more than 4 boundary lines, which gmsh does not allow. In the file I send you the ruled surface still has more than 4 boundary lines. Most fascinating, generating such a ruled surface by extrusion does not generate an error... CU andi -- Andreas Hasenkopf Phone: +49 151 11728439 Homepage: http://www.hasenkopf2000.net GPG Pub Key: http://goo.gl/4mOsM
Point(2) = {52, 0, 0, 1e+22};
Point(3) = {0, 52, 0, 1e+22};
Point(5) = {-52, 0, 0, 1e+22};
Point(6) = {0, -52, 0, 1e+22};
Point(7) = {-22, -22, 0, 1e+22};
Point(8) = {6.5, -22, 0, 1e+22};
Point(9) = {4.5, -22, 0, 1e+22};
Point(10) = {-48.5, -22, 0, 1e+22};
Point(11) = {-50.5, -22, 0, 1e+22};
Point(12) = {-22, 6.5, 0, 1e+22};
Point(13) = {-22, 4.5, 0, 1e+22};
Point(14) = {-22, -48.5, 0, 1e+22};
Point(15) = {-21, -50.5, 0, 1e+22};
Point(16) = {-23, -50.5, 0, 1e+22};
Point(17) = {-23, -52, 0, 1e+22};
Point(18) = {-21, -52, 0, 1e+22};
Point(19) = {22, 22, 0, 1e+22};
Point(21) = {22, -4.5, 0, 1e+22};
Point(22) = {22, 48.5, 0, 1e+22};
Point(27) = {50.5, 22, 0, 1e+22};
Point(28) = {48.5, 22, 0, 1e+22};
Point(29) = {-4.5, 22, 0, 1e+22};
Point(30) = {-6.5, 22, 0, 1e+22};
Point(33) = {22, 50.5, 0, 1e+22};
Point(35) = {21, -6.5, 0, 1e+22};
Point(36) = {23, -6.5, 0, 1e+22};
Point(37) = {23, -29.2, 0, 1e+22};
Point(38) = {21, -30.9, 0, 1e+22};
Point(39) = {22, 52, 0.1, 1e+22};
Point(40) = {52, 22, 0.1, 1e+22};
Point(42) = {-52, -22, -0, 1e+22};
Circle(1) = {13, 7, 9};
Circle(2) = {9, 7, 14};
Circle(3) = {14, 7, 10};
Circle(4) = {10, 7, 13};
Circle(5) = {12, 7, 8};
Circle(6) = {8, 7, 15};
Circle(7) = {16, 7, 11};
Circle(8) = {11, 7, 12};
Circle(9) = {22, 19, 28};
Circle(10) = {28, 19, 21};
Circle(11) = {21, 19, 29};
Circle(12) = {29, 19, 22};
Line(24) = {6, 18};
Line(25) = {15, 18};
Line(26) = {18, 17};
Line(27) = {17, 16};
Line(30) = {5, 3};
Circle(41) = {30, 19, 33};
Circle(42) = {33, 19, 27};
Circle(43) = {27, 19, 36};
Circle(44) = {35, 19, 30};
Line(45) = {3, 39};
Line(46) = {5, 42};
Line(47) = {6, 38};
Line(48) = {38, 37};
Line(49) = {38, 35};
Line(50) = {36, 37};
Line(51) = {37, 2};
Line(52) = {2, 40};
Circle(53) = {42, 7, 17};
Circle(54) = {39, 19, 40};
Line Loop(33) = {4, 1, 2, 3, -7, -27, -26, -25, -6, -5, -8};
Plane Surface(33) = {33};
Line Loop(37) = {12, 9, 10, 11};
Plane Surface(37) = {37};
Line Loop(38) = {4, 1, 2, 3};
Plane Surface(38) = {38};
Line Loop(56) = {12, 9, 10, 11, -44, -49, 48, -50, -43, -42, -41};
Plane Surface(56) = {56};
Line Loop(58) = {47, 49, 44, 41, 42, 43, 50, 51, 52, -54, -45, -30, 46, 53, 27,
7, 8, 5, 6, 25, -24};
Plane Surface(58) = {58};
Delete {
Surface{58};
}
Delete {
Line{45, 54, 52};
}
Delete {
Point{39, 40};
}
Point(43) = {52, 22, 0};
Point(44) = {22, 52, 0};
Line(59) = {2, 43};
Line(60) = {3, 44};
Circle(61) = {44, 19, 43};
Line Loop(62) = {8, 5, 6, 25, -24, 47, 49, 44, 41, 42, 43, 50, 51, 59, -61,
-60, -30, 46, 53, 27, 7};
Plane Surface(63) = {62};
Extrude {0, 0, 9.5} {
Surface{56, 33};
}
Extrude {0, 0, 20} {
Surface{37, 63, 38};
}
Line(329) = {230, 164};
Line(330) = {172, 200};
Line Loop(331) = {309, 310, 311, 308};
Line Loop(332) = {201, 202, 203, 204, -329, 220, 221};
Plane Surface(333) = {331, 332};
Line Loop(334) = {180, 181, 182, 179};
Line Loop(335) = {208, 209, 210, 211, 212, -330, 207};
Plane Surface(336) = {334, 335};
Extrude {0, 0, -9.5} {
Surface{336, 333};
}
Point(338) = {-22, -46.5, 0};
Point(339) = {-22, 2.5, 0};
Point(340) = {-46.5, -22, 0};
Point(341) = {2.5, -22, 0};
Point(342) = {-2.5, 22, 0};
Point(343) = {46.5, 22, 0};
Point(344) = {22, -2.5, 0};
Point(345) = {22, 46.5, 0};
Circle(451) = {345, 19, 342};
Circle(452) = {342, 19, 344};
Circle(453) = {344, 19, 343};
Circle(454) = {343, 19, 345};
Circle(455) = {339, 7, 340};
Circle(456) = {340, 7, 338};
Circle(457) = {338, 7, 341};
Circle(458) = {341, 7, 339};
Line Loop(459) = {451, 452, 453, 454};
Plane Surface(460) = {459};
Line Loop(461) = {455, 456, 457, 458};
Plane Surface(462) = {461};
Extrude {0, 0, -20} {
Surface{460, 462};
}
Point(372) = {-21, -52, -20};
Point(373) = {-23, -52, -20};
Point(374) = {0, -52, -20};
Point(375) = {-52, -22, -20};
Point(376) = {-52, 0, -20};
Point(377) = {52, 0, -20};
Point(378) = {0, 52, -20};
Point(379) = {22, 52, -20};
Point(380) = {52, 22, -20};
Point(381) = {21, -30.9, -20};
Point(382) = {23, -29.2, -20};
Line(507) = {379, 378};
Line(508) = {378, 376};
Line(509) = {376, 375};
Line(510) = {373, 372};
Line(511) = {372, 374};
Line(512) = {374, 381};
Line(513) = {381, 382};
Line(514) = {382, 377};
Line(515) = {377, 380};
Circle(516) = {379, 347, 380};
Circle(517) = {375, 360, 373};
Line Loop(518) = {464, 465, 466, 467};
Line Loop(519) = {489, 486, 487, 488};
Line Loop(520) = {508, 509, 517, 510, 511, 512, 513, 514, 515, -516, 507};
Plane Surface(521) = {518, 519, 520};
Extrude {0, 0, 19} {
Surface{521};
}
Point(463) = {22, 46.5, 20};
Point(464) = {22, -2.5, 20};
Point(465) = {-2.5, 22, 20};
Point(466) = {46.5, 22, 20};
Point(467) = {-46.5, -22, 20};
Point(468) = {2.5, -22, 20};
Point(469) = {-22, 2.5, 20};
Point(470) = {-22, -46.5, 20};
Circle(619) = {468, 149, 470};
Circle(620) = {470, 149, 467};
Circle(621) = {467, 149, 469};
Circle(622) = {469, 149, 468};
Circle(623) = {466, 136, 464};
Circle(624) = {464, 136, 465};
Circle(625) = {465, 136, 463};
Circle(626) = {463, 136, 466};
Line Loop(627) = {626, 623, 624, 625};
Plane Surface(628) = {627};
Line Loop(629) = {622, 619, 620, 621};
Plane Surface(630) = {629};
Extrude {0, 0, 20} {
Surface{628, 630};
}
Point(497) = {-21, -52, 40};
Point(498) = {-23, -52, 40};
Point(499) = {0, -52, 40};
Point(500) = {-52, -22, 40};
Point(501) = {-52, 0, 40};
Point(502) = {52, 0, 40};
Point(503) = {52, 22, 40};
Point(504) = {22, 52, 40};
Point(505) = {0, 52, 40};
Point(506) = {21, -30.9, 40};
Point(507) = {23, -29.2, 40};
Line(675) = {503, 502};
Line(676) = {502, 507};
Line(677) = {507, 506};
Line(678) = {506, 499};
Line(679) = {499, 497};
Line(680) = {497, 498};
Line(681) = {500, 501};
Line(682) = {501, 505};
Line(683) = {505, 504};
Circle(684) = {503, 472, 504};
Circle(685) = {498, 485, 500};
Line Loop(686) = {633, 634, 635, 632};
Line Loop(687) = {655, 656, 657, 654};
Line Loop(688) = {678, 679, 680, 685, 681, 682, 683, -684, 675, 676, 677};
Plane Surface(689) = {686, 687, 688};
Extrude {0, 0, -19} {
Surface{689};
}
#!/usr/bin/env python
import pylab
import fipy
import scipy.integrate as si
# 3D mesh
mesh = fipy.GmshImporter3D('twin_ap_3D_X4.msh')
X,Y,Z = mesh.getFaceCenters()
z0,r,b = 10.,27.5,1.
zb,zt = -0.5,20.5
x0,y0 = 22,22
x1,y1 = -22,-22
potential = fipy.CellVariable(mesh=mesh, name='potential', value=0.)
permittivity = 1.
potential.equation = (fipy.DiffusionTerm(coeff = permittivity) == 0.)
allfaces = mesh.getExteriorFaces()
# 3D faces
ring0 = allfaces & ( (( ((X-x0)**2+(Y-y0)**2 > (r-b-.5)**2) & ((X-x0)**2+(Y-y0)**2 < (r+b+.5)**2) ) | ((abs(X-x0)<b+.6) & (Y<y0-r-b))) & (abs(Z-z0)<.6) )
ring1 = allfaces & ( (( ((X-x1)**2+(Y-y1)**2 > (r-b-.5)**2) & ((X-x1)**2+(Y-y1)**2 < (r+b+.5)**2) ) | ((abs(X-x1)<b+.6) & (Y<y1-r-b))) & (abs(Z-z0)<.6) )
bottom_faces = allfaces & ( (((X-x0)**2+(Y-y0)**2 > (r-b-2.1)**2) ) & ( ((X-x1)**2+(Y-y1)**2 > (r-b-2.1)**2) ) & (abs(Z-zb)<1.) )
top_faces = allfaces & ( (((X-x0)**2+(Y-y0)**2 > (r-b-2.1)**2) ) & ( ((X-x1)**2+(Y-y1)**2 > (r-b-2.1)**2) ) & (abs(Z-zt)<1.) )
the_top = allfaces & (Z>39.99)
bcs = (
fipy.FixedValue(value=5.,faces=ring0),
fipy.FixedValue(value=0.,faces=ring1),
fipy.FixedValue(value=1.,faces=bottom_faces),
fipy.FixedValue(value=0.,faces=top_faces),
# fipy.FixedValue(value=0.,faces=mesh.getFacesLeft()),
# fipy.FixedValue(value=0.,faces=mesh.getFacesRight()),
# fipy.FixedValue(value=-1.,faces=mesh.getFacesTop()),
fipy.FixedValue(value=-1.,faces=the_top),
)
potential.equation.solve(var=potential, boundaryConditions=bcs)
print ">>> Solved <<<"
# VTKViewer for 3D mesh
viewer = fipy.viewers.vtkViewer.VTKViewer(vars=(potential,))
viewer.plot("/tmp/test.vtk")
print ">>> FINISHED <<<"
signature.asc
Description: OpenPGP digital signature
_______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
