Add support for facet orientation in FacetBasis #865
Conversation
|
Hmm. I think the only times that i've used an internal curve (in 2-D) were for lineal sources, like point sources distributed along a curve. For these the sign of the normal doesn't matter; it's just a linear form like I suppose cases where the sign of the normal is required might arise, though I can't think of any now. Maybe if the internal interface is specified as a sequence of facets (in two dimensions) then a consistent orientation could be inferred from that but obviously that's impossible in three dimensions. |
By default (using |
|
I'm referring to your comment in #854:
|
Right. For n = 2, FreeFEM includes commands for defining two-dimensional geometry and meshing it with triangles. First one defines curves parametrically, like (x (t), y (t) for a ≤ t ≤ b and then these curves can either be closed or chained together to bound a planar region. These curves can be tagged with an |
For a lineal source? (Or a superficial source in three dimensions?) No, I don't think it's an issue, because the normal isn't referred to in the linear form, so I don't think it matters if the sign of the normal is arbitrary. |
Alright, so there is some property in msh files that tells the orientation and we need to understand. |
|
Well, yes, probably, for the case of an internal surface in three dimensions, but that's not terribly common so i wouldn't think it urgent. I don't think i've ever actually used that in 3-D myself. Two dimensions yes, but that's easier because it's possible to specify the facets in sequence. I think the details of how the orientation is specified by the user can be left open for now? On the understanding that users might want to specify it in different ways. A really common way is that the defining facets are specified not as facets of cells but by their points, with the order of the points significant. I tripped over this way back in #72. It had not occurred to me that meshes might not imply the orientation of cells. Anyway, this is how Gmsh does it on MSH: facets are (n−1)-dimensional cells supported on their own points with the sequence implying orientation. |
|
Yeah, that might be a reasonable thing to support: optionally inferring the orientation of the "boundaries" if reading from meshio. I think MSH is not alone in linking the sequence of points with orientation. Typically it's a right-hand rule for trianglar facets. For line-facets in two dimensions, it's the same: the line is assumed to be described with its normal pointing right; as for an outward normal when a boundary is described counterclockwise. |
|
Our description of facets is more implicit though:
Why I wrote "yet"? I've been thinking about the integration of bilinear forms over 2D surfaces embedded in 3D. There you would need an explicit description of the surface. I've done some work towards supporting this but I'd like to do it properly before releasing. What do I mean by "properly"? There are examples where you have one unknown on, say, a tetrahedral mesh and another unknown on it's boundary. An interesting example is elasticity equations (displacement) on a tetrahedral mesh and Reynolds equation (pressure) on the boundary. Now while you can sort of hack this together already using However, this seems to be somewhat nontrivial because of how I have several ideas and are just trying to think about the best course of action. However, this is unrelated to what we are doing here so I'll stop the rant now. |
|
Thank you; super interesting. Is this a bit like #514 ? There I wanted what I called there a FacetMesh, to be extracted from a MeshQuad (but conceptually any skfem.Mesh of There are also applications involving surfaces in space not attached to a volume mesh, e.g. minimal or constant-mean-curvature surfaces which cannot be expressed as z = z (x, y) as in ex10); I have solved such problems using Surface Evolver. These are deforming curved surface meshes. I suppose curves in the plane and space are also possible, e.g. two-dimensional or axially symmetric reductions of the former. |
Yes, I think that would be one of the use cases. Later this year I'm traveling abroad for a research visit where the plan is to study "shell models" using differential geometry. So solving equations on embedded surfaces will definitely be a short term priority for me. |
|
If I'm understanding what @gdmcbain is saying, sometimes he wants flux across a subset of facets that don't necessarily form the boundary of a subdomain. I think maybe class AbstractBasis: ...
class CellBasis(AbstractBasis): ...
class FacetBasis(AbstractBasis): ...
class InteriorFacetBasis(FacetBasis): ...
class BoundaryFacetBasis(FacetBasis): ...
class InteriorBoundaryFacetBasis(BoundaryFacetBasis): ...And I guess when I write it this way, I can understand how #865 I wish we could get rid of class BoundaryFacetBasis(FacetBasis):
def __init__(self, ..., elements, facets, outward_normal, outward_trace):
...If |
|
Part of what I'm thinking is that in software terms none of the subclasses of have any methods other than If my understanding there is right, then that refactoring can be considered and accepted or rejected without too much controversy. The important thing then, either way, is how to build the subsets of facets along with their orientation, either way, whether that's in the There's a lot of lines of code in this PR, but if scikit-fem/skfem/assembly/basis/boundary_facet_basis.py Lines 105 to 106 in 85bea6c is the core of it and the rest is just along for the ride, then yeah I think this is what I'm talking about and I think it's what's wanted, along with possible siblings to be added to the library or in user-level code as required. By 'unpythonic' do you mean the name? I'm not terribly fond of Yeah, rereading your comment, I think maybe we're converging? I think the name |
|
I think we're converging too. Let me just add this analogue to clarify what I'm saying. If we made a geometry module, they might all be polygons. And you could argue that we should just have one But I argue that by subclassing Polygon into Square and Triangle, we get to make more, simpler constructors (the unix mentality). We also get a convenient place to do extra parameter checks, raise informative exceptions, place extra documentation specific to those shapes. This is purely a semantic/syntactic computer sciency argument, and it can be taken to far. 100's of classes for every named Polygon would become unwieldy, hard to remember the right name, and most of that code would never get run in userland, ever. But for AbstractBasis, even with the feature creep and the extensions we're proposing now, its still less than 10 classes. Quite manageable, and they're all really useful. Keeping with this analogy, by unpythonic, I meant that I dislike class Polygon:
@classmethod
def init_square(...):...
def init_triangle(...):...
def __init__(self,...):...Nothing stops from writing the class this way, but it doesn't feel like the spirit of python or object oriented to me. It's kind of like having the naked A last point about it, is that subclassing for Square and Triangle is going to be a familiar construction for many python programmers. classmethods for init_square, etc, are not going to be familiar, and consequently they will steepen the learning curve for new adoptors. I think there's a lot of value in staying as close to familiar idioms as possible. To answer your questions:
Whether it becomes unnecessary or not, the ability to tweak the methods of AbstractBasis before calling super() is one of the top reasons why I prefer subclassing over fancy constructors. It usually makes simpler, understandable code with less if statements.
And particular exceptions and particular documentation, and occasionally specialty functionality / attributes. Example: Square might implement "diagonal_length()" which might not make sense for other things. I don't have an example at hand right now for skfem, but I feel sure one will turn up. (is_boundary_closed comes to mind, but meh.)
Same thing, as a general rule, I'm opposed to large functions with lots of if statements. Sometimes this rule is hard to keep (get_dofs, for example) but in this case particularly, it is one of the main use cases for subclassing.
No, but by now you know what i meant :)
Where possible, I prefer simpler functions with less if statements. Optional keyword arguments gets very iffy.
Computationally, I agree. Conceptually, no, its very relevant because as a user you're trying to do something, and when it doesn't work, you want messages like "Exception: can't create cell basis from a facet boundary that isn't closed"... not messages like "Exception: None type has no method is_closed" |
|
@gdmcbain thoughts? are we converging? class AbstractBasis:
def __init__(self, mesh, elem, etc):
...
class CellBasis(AbstractBasis):
def __init__(self, mesh, elem, etc, elements=None):
"""Define a basis over a set of cells.
elements
an arbitrary set of elements in the mesh
For `elements`, None is interpreted as "the set of all elements"
"""
...
class FacetBasis(AbstractBasis):
def __init__(self, mesh, elem, etc, facets=None):
"""Define a basis over a set of facets. Use subset of mesh if facets is not None.
facets
an arbitrary set of facets in the mesh.
For `facets`, None is interpreted as "the set of all facets"
"""
...
class InteriorFacetBasis(FacetBasis):
def __init__(self, mesh, elem, etc, elements=None, facets=None, combine="intersect"):
"""Define a basis over a set of facets. Excludes facets that bound a single cell of the set specified by `elements`.
elements
specifies the set of facets interior to these elements.
facets
an arbitrary set of facets in the mesh.
combine
specify how to combine the sets specified by `elements` and `facets`
"union" - union of the two sets
"intersect" - intersection of the two sets
For `facets` and `elements`, None is interpreted as "the set of all facets"
"""
...
class BoundaryFacetBasis(FacetBasis):
def __init__(self, mesh, elem, etc, elements=None, facets=None, combine="intersect"):
"""Define a basis over a set of facets. Excludes facets that bound two cells of the set specified by `elements`.
elements
specifies the set of facets that bound these elements.
facets
an arbitrary set of facets in the mesh.
combine
specify how to combine the sets specified by `elements` and `facets`
"union" - union of the two sets
"intersect " - intersection of the two sets
For `facets` and `elements`, None is interpreted as "the set of all facets"
"""
... |
|
No, the idea that a class's having multiple constructor functions is unpythonic is new to me; e.g. for
Afterwards one just has an |
No. In Unix, every command has the same type of output, line-separated text, and every command returns the same type value, an int (0 for no problem). As Alan Perlis put it in the ninth of his epigrams on programming:
I suspect the success of the ecosystem growing up around NumPy is in large part due to the use of
There was an ugly wobble with
For the polygonal example, I do quite a bit of work with polygons (e.g. integrated circuit layout) and usually it is just a single class for n-gons for n > 2 rather than a denumerably infinite sequence of subclasses But yes, resoundingly, to ‘more, simpler constructors’; this is indeed the Unix philosophy, and Perlis's, and, I think, is exactly what we want. Whether they're classmethods like this PR or module-functions, I don't think matters much; whatever makes for simpler implementation in each case. The classmethod here looks fine. |
|
These are all really good points, and I threw "unpythonic" around too flippantly. And in particular, your invocation of numpys maxtrix is a powerful one. I too am grateful that matrix type died. I think I overreached in this PR, because your counterexamples are compelling. But I'm still quite interested in helping make skfem easier to use if I can. Are you proposing collapsing |
Yes, I think so. I need to study
Well the obvious difference initially is that facets are (n − 1)-dimensional and so have a normal and a side (for the traces). But lately @kinnala mentioned the very exciting extension to formulating problems on what are now facet-bases. There was a little bit of this in #514 (just L2-projections on rectilinear facet-bases) but if it will mean that general meshes have nontrivial metric tensors then I suppose maybe there will be less difference between an n-dimensional cell-basis and an (n − 1)-dimensional facet-basis? That is, both might equally well be curved with respect to some still higher-dimensional space that they're embedded in. I guess the main thing about a facet-basis though is that is related to some cell-basis of one-higher dimension. The facets which are the cells of the facet-basis are facets of cells of the cell-basis. And they're linked by trace operators: if there's a function defined on the cell-basis then it has values (traces) implied on the facet-basis. Also it has the normal component of the gradient, which is another trace operator. There are probably higher ones too, but these are the two most met in applications governed by second-order partial differential operators. Thus we do have |
maybe irrelevant to this PR, but I'm extremely interested in this topic. might be a good general discussion to start. i'm trying to embed skfem in an impedance tomography problem, and assumptions about the electrode interface is pretty significant to the found solutions.
e.g. how the latitudes are curved on the surface of a sphere, which is itself curved in 3d? I can't seem to stop thinking about (what I think of as) pathological meshes... e.g. mobius strip... and how these concepts of orientation apply. It's very non-trivial for me to try to visualize what "consistent normal" might look like this n-dimensional context.
I still want to pull one thing from my "unpythonic" rant... I'm specifically thinking about "should basis.zeros() return an ndarray or a subclass"...
I said "no" before, but i want to amend that to "sometimes". the But having read all this the past couple of days, I'm in your camp now that CellBasis and FacetBasis seem like enough. Whether AbstractBasis is enough has gotten abstract enough that I can't really visualize it :) |
Great. I launched #867 ‘interfacial coupling and discontinuities’. |
Yes, exactly. If the curve or surface can be described parametrically in coordinates in the embedding space, e.g. x, y, and z as functions of latitude and longitude for a sphere, or functions just of longitude for a parallel, then there are formulæ for calculating the curvature, normalls, binormals, torsion, &c., from these. These may be found in treatises on differential geometry, particularly Riemannian geometry. The key concept is the metric tensor which in simple terms gives the net infinitesimal distance implied by an infinitesimal step in each of the local coordinates. There are also formulæ for the case in which the surface is given by a level-set of a function in space; here the normal is parallel to the gradient of that function and the mean curvature is given by the divergence of the gradient. It's also possible to start with a metric tensor and not think of the working space as being embedded in a higher-dimensional space; this might arise in relativity (and so in some formulations of electromagnetism?) but not so much in my sublunary applications so far. |
Right. A Möbius strip is nonorientable as a surface. It could be triangulated easily enough but the resulting mesh would not be orientable either. I think most meshers would baulk; e.g. Gmsh. By far the most dealings that I have had with inconsistently oriented manifolds apart from #72 and #821 have been in exceptions raised by Gmsh because I have tried to define a volume by its bounding surface which I have stitched together from multiple patches with inadvertently inconsistent orientations but, distinguishing between oriented and orientable, I have never had to deal with nonorientable manifolds and am not sure how they might arise or even what kind of boundary value problems might reasonably be posed on them. It'd be easy enough to parameterize the Möbius strip as embedded in ℝ³ and then try to calculate the normal using the standard formulæ. I haven't tried that. I assume it's two-valued? There's probably a square-root in there somewhere?
Yeah, it's best to fall back on the standard formulæ from differential geometry which visualization fails or deceives, then build up the intuitiion again from there. The introductions to this field that I found particularly useful and intuitive were:
but there might well be more modern ones that are good too. Very likely there are even better ones relating to electromagnetism; see https://en.wikipedia.org/wiki/Differential_form#Applications_in_physics. |
Hmm, I see. Yeah, I'm used to just remembering which scikit-fem/docs/examples/ex29.py Lines 87 to 90 in b24b078 So I think of |
I'd like to understand |
|
Yeah nah, a global sign attribute seems like a cumbersome effort to avoid minus-signs in formulæ which is impossible: move a term from one side of an equation to the other and it changes sign, e.g. Please show a plausible case in which this would result in simpler code. Conditionals inside forms öook reslly bad. |
That is done using FacetBasis(..., side=0/1) which does not affect the normal. |
numpy is this way too, I think. I write It is possible in numpy to change the dtype of an array you already have. Might I construct new facet basises this way? Similar to b1 = FacetBasis(...)
b2 = b1.with_flipped(normals=False, traces=True)
I hadn't thought of it this way, but I actually like the metaphor a lot. I'm not sure I'm giving good use cases, but it does seem like the facet basises are "signed", and if so, I think you should be able to get at the sign? Maybe I'm not proposing the right implementation. |
Oh... maybe I'm misunderstanding the new syntax then. |
Is this because the design goal of a forms/functional is one line of code that looks very similar to what's written on blackboard? |
I haven't pulled your latest commits, but I did move from #851 to #865. Currently I'm using f_p2 = skfem.FacetBasis(mesh, skfem.ElementTriP2(), facets=mesh.facets_around('inner'), flip_traces=True)This is giving me outside traces and outside pointing normals. a = FacetBasis(m, e, facets=m.facets_around('subdomain'), flip=False, side=0) # default
b = FacetBasis(m, e, facets=m.facets_around('subdomain'), flip=False, side=1)
c = FacetBasis(m, e, facets=m.facets_around('subdomain'), flip=True, side=0)
d = FacetBasis(m, e, facets=m.facets_around('subdomain'), flip=True, side=1)Am I just doing things all wrong to be worrying so much about which trace and which way the normal? It is currently very important to keep track of in my code. In most cases I care about, the traces aren't equal. Often (but not always) I have a DC boundary condition on one side, and I never want that trace. A lot of my bugs have come down to which way the normals are and which trace I'm using. edit: Almost everything I'm doing with skfem involves working with net fluxes into and out of specific regions in the domain. |
Fixed. |
Is it possible to use Gauss divergence theorem and turn it into a domain integral? |
|
This PR is hopefully ready soon. I think I'll need to leave the save/load later because I'm running out of time. |
No and yes. The trace is really important because it can have completely different values on the two sides (e.g. ex26); it does need to be worried about and selected carefully. Yes for the normal though; as long as it's consistent, it's just a matter of global sign and this cannot be fixed with an attribute. It might be reasonable to define the unary negation of an oriented facet array as the same array with all the orientations flipped. This is consistent with the mathematical idea of a contour integral in the plane with a counterclockwise or clockwise contour; similarly for surface integrals in space. |
Yeah, I wonder whether we were a bit too clever in #681 ...or do we just need a little bit more of such cleverness here. I have not successfully run the failing encode decide tests in the mental Python interpreter or gotten to a terminal so I'm not yet clear on the issue. I hope to get to this next week but I would think it reasonable to proceed with the rest of the PR as serialization is fairly specialized. |
Almost certainly, and this will be more robust and accurate as demonstrated in our ANZIAM Journal paper and ex. 13. |
Discussion: #873 |
|
Oops! I almost missed these comments. Does github preserve these discussions on deleted branches indefinitely? There's a lot in this one I want to try to migrate over to more formal documentation when I have time.
I think this is along the lines of what I was trying to describe. And that idea of being able to In any case, @kinnala set me straight on the syntax, and so now all my facet basises are just born correct and I don't foresee a reason to want to be flipping the normals around or ever testing which way they go.
I'll look into these, thanks for the reference. |
Is this the relevant portion? This is basically what I am doing, except my conductivity is not uniform. I thought you meant to change from an integral over facets to an integral over cells. scikit-fem/docs/examples/ex13.py Lines 55 to 64 in 33a2513 |
|
If you apply Gauss divergence on that integral around a closed loop you can replace it by u.T @ A @ u where A is the matrix (sigma grad(u), grad(v)) integrated over the subdomain. But I think this requires u is not constant inside the subdomain? |
|
No, that's a boundary integral; the bit in ex13 that's being recommended is scikit-fem/docs/examples/ex13.py Line 51 in 33a2513 It doesn't immediately look like a 'domain integral', more like a quadratic form, but this is the third of the ''Three ways to compute multiport inertance'; the boundary flux integral is the first and the connection is made explicit by the second. |
Oh yeah you'd want to evaluate the (sub)domain integral over the complementary subdomain in which the potential actually varies between the source and sink values. |
The nonuniformity shouldn't matter, even anisotropy, so long as the conductivity is everywhere a positive definite symmetric tensor bounded above a positive constant. It might be nice to generalize the derivation in the ANZIAM J. to this case; I think that should be straightforward. Insofar as I understand the application it's going to conceptually involve calculating the conductance matrix for a bunch of electrodes ('ports') embedded in a conducting solid and there'll be one such matrix for each geometry, set of electrode positions, and conductivity fields. |
|
What the conductance matrix is is that if one charges the ith electrode to unit potential and earth the rest, the ijth coefficient is minus (can't avoid those minus signs!) the current out of the jth electrode. |
|
Further generalization: if one uses AC instead of DC the conductivity becomes complex and the symmetric conductance matrix becomes a Hermitian admittance matrix. It should still be positive though if the medium is passive. |
No. A |
|
I don't think it's really worth defining the unary negation of an |
Fair. After I wrote that comment it did occur to me that even just a single complex number has unary negation but there's no good concept of the "sign" of a complex number. So clearly not all things that can be negated are signed.
I agree completely. After the dust settled I think FacetBasis and OrientedBoundary came out making a lot of sense. All that sign/negation stuff was born out of some confusion mid-PR. |
This PR now acknowledges the existence of oriented boundaries/interfaces (i.e. sets of facets).
We can load Gmsh file with an oriented interface in the middle:
Here
OrientedBoundaryis a subclass ofndarraywith the additional attributeoriwhich is either 0 or 1 for each facet.We can visualize the orientation:
The orientation is taken into account in
FacetBasis, e.g.,Obviously there can be multiple oriented facet sets:
There is a facility for creating oriented facet sets around subdomains in
Mesh.facets_around:We can also now orient facet sets constructed with
Mesh.facets_satisfying:Fixes #821 and fixes #870.
Missing work: boundary orientations are invalidated by save/load cycle and by refine.