Towards XFEM - Cut Cell Quadrature

[benkirk]

No rush on this, but this is 40 commits already and no one has seen it, so I'm putting this out there for some initial feedback.

The API is basically demonstrated in the new examples/subdomains/subdomains_ex3 where we seek to integrate a discontinuous function not aligned with the mesh.

The approach is as follows: QComposite defines a composite quadrature rule based on an underlying Gauss quadrature rule.

For uncut cells, this reverts to a standard Gauss rule, for cells cut by an interface, the cell is decomposed into simplices via the new ElemCutter class, whose current implementation uses triangle/tetgen, and agglomerates the simplex quadrature points/weights into a "composite' rule.

The current implementation requires being able to provide COMM_SELF to the "working" meshes created inside the ElemCutter class.

The ElemCutter class right now takes vertex values of a signed distance function to locate the cut, but this is ripe with opportunity for extension.

[benkirk]

@jwpeterson - have you looked at all at QComposite in here?

As the father of pretty much all our quadrature, I'd be especially interested in your thoughts...

[jwpeterson]

On Apr 13, 2013, at 7:44 AM, "Benjamin S. Kirk" notifications@github.com wrote:

@jwpeterson - have you looked at all at QComposite in here?

As the father of pretty much all our quadrature, I'd be especially interested in your thoughts...

Sure thing... In a broad sense, I've heard it's better to use low-order quadratures with more intervals when dealing with discontinuous functions...

[roystgnr]

QComposite uses the correct subcell decomposition when combining its
quadrature rules, right? So it should still be effective to use
higher order rules.

In the general case of a discontinuity of unknown location, though,
higher-order rules are indeed pointless and the general rule is to
bound the discontinuity position in as tight an interval as possible.

[benkirk]

On Apr 13, 2013, at 11:02 AM, "roystgnr" notifications@github.com wrote:

QComposite uses the correct subcell decomposition when combining its
quadrature rules, right? So it should still be effective to use
higher order rules.

That's right. When a distance function is present, it resolves the discontinuity exactly and used whatever rule you requested on the subcells whose union defines the two parts of the cut element.

There is an important partial specialization missing right now - for various netwon-cotes rules the subcells will produce duplicate points which are harmless but should be detected and combined for matrix assembly efficiency.

In the general case of a discontinuity of unknown location, though,
higher-order rules are indeed pointless and the general rule is to
bound the discontinuity position in as tight an interval as possible.

My thoughts here would be to use the cut cell stuff when the interface location is known, and perhaps Monte Carlo or something when it is not?

-Ben

[buildhive]

libMesh - C++ Finite Element Library » libmesh #490 SUCCESS
This pull request looks good
(what's this?)

[jwpeterson]

So if I understand correctly, the main ideas of QComposite are:
1.) QBase::init() now becomes a virtual, overloaded method, and takes a const std::vector &vertex_distance_func.
2.) QComposite is a template class which inherits from its template argument (a version of the curiously recurring template pattern?)

I guess #2 is a reasonable approach vs. just holding a QBase* in QComposite and forwarding relevant calls to it?

One downside of the template approach is that you have to do explicit instantiations for QComposite et al. and add all the 'using QSubCell::xyz' stuff.

Finally, would you ever be in a situation where you'd need QComposite because then QComposite would be derived from itself ;-)

Some other comments:

I've found that tetgen is generally a pretty bad piece of software... From your comments, it looks like it's generating zero volume elements for you? Have you tried upgrading our version? Looks like we have 1.4.1 but there may be a 1.4.3 available. I'm also curious if there's anything better out there...

Just curious, what's the idea behind the numbering of our recent quadrature rule enumerations: QCLOUGH=21, QCOMPOSITE=31?

[jwpeterson]

Oops, markdown. That should have read "...explicit instantiations for QComposite<QTrap> et al."

[jwpeterson]

And "...would you ever be in a situation where you'd need QComposite<QComposite> because then QComposite would be derived from itself ;-)"

Totally ruining the joke...

[benkirk]

I've found that tetgen is generally a pretty bad piece of software... From your comments, it looks like it's generating zero volume elements for you? Have you tried upgrading our version? Looks like we have 1.4.1 but there may be a 1.4.3 available. I'm also curious if there's anything better out there...

I'll see if it is worth an upgrade. Notably, when allowing tetgen to generate 'quality' meshes it does not produce the zero volume cells, but it does produce way more than I need. Since for this application the interpolation properties of the subtets are irrelevant (all they are used for is mapping quadrature points) its not such a big deal.

Longer term, though, if this becomes a commonly used funcitonality we could perhaps reimplement this to use QHull or something else.

Just curious, what's the idea behind the numbering of our recent quadrature rule enumerations: QCLOUGH=21, QCOMPOSITE=31?

@roystgnr - any comments on the 21? I was just kinda following suit, for no good reason other than maybe someone might want to mod it with 30 - although I'd highly discourage that kind of thing!

[roystgnr]

Whoever added our Infinite element types to the FEFamily enum added a
bit of space in there for expansion, so I did the same, and then I got
carried away and did the same for quadrature.

We could actually reorganize the quadrature enums if we wanted, I
think, since those numbers are just internal to each run. I think we
write FEFamily out to file, though, so there backwards compatibility
leaves us stuck with whatever initial numbering each entry gets.

[buildhive]

libMesh - C++ Finite Element Library » libmesh #528 SUCCESS
This pull request looks good
(what's this?)

[benkirk]

@jwpeterson wrote

I've found that tetgen is generally a pretty bad piece of software... From your comments, it looks like it's generating zero volume elements for you? Have you tried upgrading our version? Looks like we have 1.4.1 but there may be a 1.4.3 available. I'm also curious if there's anything better out there...

Yes, I've got a similar impression. I'm looking into Qhull as an alternative: http://www.qhull.org

I think it could easily do everything I'm asking of Tetgen in this particular application - namely, return the convex hull of a small number of points in 3D.

-Ben

[permcody]

On Wed, May 8, 2013 at 8:59 AM, Benjamin S. Kirk
notifications@github.comwrote:

@jwpeterson https://github.com/jwpeterson wrote

I've found that tetgen is generally a pretty bad piece of software... From
your comments, it looks like it's generating zero volume elements for you?
Have you tried upgrading our version? Looks like we have 1.4.1 but there
may be a 1.4.3 available. I'm also curious if there's anything better out
there...

Yes, I've got a similar impression. I'm looking into Qhull as an
alternative: http://www.qhull.org

I think it could easily do everything I'm asking of Tetgen in this
particular application - namely, return the convex hull of a small number
of points in 3D.

-Ben

That's interesting! At one point I was debating on using a 3D hull for one
of the problems I'm working on, but I decided not to go that direction
because I didn't have a readily available package to do so, and I believed
that using a simple bounding sphere would be more effective for my problem.
In my case I need to compute intersections of "bodies" but it's a fuzzy
comparison so spheres are working nicely for my problem. I like the idea
of having this available if I want to try something else later though.

Cody

—
Reply to this email directly or view it on GitHubhttps://github.com//pull/84#issuecomment-17611381
.

[cato-]

Is it also possible to use this integration with other cutting geometries than just linear level-set-functions?
A similar integration is also needed in mortar contact codes where you need to integrate over the exact contact interface. This boundary of this interface might also go through an element and is in general not just a line or plane like in the XFEM case.

[jwpeterson]

This PR needs some freshening up. I'll reopen a new PR for the rebased branch, and refer back to this one.