On Dec 21, 2011, at 7:27 PM, James Snyder wrote: > On Wed, Dec 21, 2011 at 5:22 PM, James Snyder <[email protected]> wrote:
>> I expect that the current measurements will be quite sensitive to the >> cell volumes, but the discontinuity mentioned seems to persist at 0.39 >> vs 0.4 even with varying mesh resolution. Am I doing something >> foolish or is there something simple I've missed here? Is the issue >> related to something happening in Gmsh or FiPy? I see that the cell >> count sometimes drops a bit at the transition, but if I include 0.38 >> and 0.42 the discrepancy doesn't seem proportional to the drop. > > Scratch this first question. I figured out what was happening. My > boundary selection parameters were including the edge of the outer > volume once I dropped below 0.40. I'd still like suggestions or > responses for the remaining questions though, if possible. We're glad you were able to sort this out. >> If I want current density to be zero at certain boundaries normal to >> that surface (insulators), what is the correct way to do that? I see >> the new constrain approach, but I'm just setting the faceGrad value >> vector to all zeros: "potential.faceGrad.constrain([[0], [0], >> [0]],outer_edges)" but this isn't quite right. Zero flux is the default boundary condition in the finite volume method, so you shouldn't need to do anything to achieve that. If you want to constrain the normal current density to some non-zero value, then see the FAQ "How do I apply a fixed flux boundary condition?" Basically, you add a source that consists of the divergence of the flux you want: >>> eqn = (TransientTerm() + ConvectionTerm(convCoeff) ... == DiffusionCoeff(diffCoeff) ... + (mesh.exteriorFaces * fluxMagnitude * mesh._faceNormals).divergence) >> I tried using physical surface labels in 3D to label faces, and while >> the labels came through, none none of the faces are actually tagged in >> mesh.physicalFaces. I'm not that interested in tagging values, but >> since one way to have "internal boundaries" is to carve out a section >> of a mesh, then select specific faces of that internal hole, it would >> be nice to use the physical labels. This is just speculation, but the documentation for "Physical Surface" at http://geuz.org/gmsh/doc/texinfo/gmsh.html#Surfaces says "the expression-list on the right hand side should contain the identification numbers of all the elementary surfaces that need to be grouped inside the physical surface." It appears to me that a Surface Loop is not an "elementary surface". If I replace your definitions with: Physical Surface("Pod") = {Ecf1[], Ecf2[], Ecf3[], Ecf4[],Bod1, Bod2, Bod3, Bod4, Ecb4[], Ecb3[], Ecb2[], Ecb1[] }; and Physical Surface("Tank") = {1002, exbox[1], 1033, exbox[6], exbox[11], exbox[16] }; then mesh.physicalFaces seem to be defined with non-empty masks. mesh.physicalCells["Water"] was already defined in your script, but mesh.physicalCells["Pod"] was empty. I noted the difference was that "Water" was defined by Volume(10014) = {10006, 10007}; Physical Volume("Water") = 10014; but "Pod" was defined by Surface Loop(10006) = {Ecf1[], Ecf2[], Ecf3[], Ecf4[],Bod1, Bod2, Bod3, Bod4, Ecb4[], Ecb3[], Ecb2[], Ecb1[] }; Physical Volume("Pod") = 10006; Changing this to Surface Loop(10006) = {Ecf1[], Ecf2[], Ecf3[], Ecf4[],Bod1, Bod2, Bod3, Bod4, Ecb4[], Ecb3[], Ecb2[], Ecb1[] }; Volume(10012) = {10006}; Physical Volume("Pod") = 10012; seems to resolve that issue, too. For the record, I don't have a clear sense of when you can reuse Gmsh identification numbers and when you can't. I think they only need to be unique in each geometry category (surfaces vs. volumes). >> Is there an easy way to get values for a cross-sectional plane of cell >> values on an unstructured mesh like the rectangular prism used in this >> example? I see that I can do quite a bit of manipulation in Mayavi, >> but it would be nice to have access to something like this without an >> extensive custom routine or jumping out to a viewer/postprocessor. There are two approaches I would take: One is to use the __call__() method of the field you are interested in to obtain interpolated values. The following constructs a 2D mesh that makes an XZ slice down the middle, then uses the coordinates of that 2D mesh to extract the the values of the 3D field at Y=0. sliceMesh = (Grid2D(Lx=tankwidth, nx=100, Ly=tankheight, ny=100) + [[-tankwidth/2.], [-tankheight/2.]]) sliceCenters = sliceMesh.cellCenters.value sliceCenters = numerix.insert(sliceCenters, 1, 0., axis=0) slice = CellVariable(mesh=sliceMesh, value=potential(sliceCenters)) The other is to extract the values of a slab within some distance of the position you are interested in: x, y, z = mesh.cellCenters slab = potential[(z < HEIGHT + THICKNESS / 2.) & (z > HEIGHT - THICKNESS / 2.)] The first is probably more useful for display purposes, the second is much faster. >> Also, Is there any likelihood or possibility that one might be able to >> do parallel computation on unstructured meshes using Trillinos in the >> future? I see that currently the parallel computation is only >> supported for gridded meshes. parallel computation on Gmsh-generated meshes is available in the next (imminent) release. _______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy]
