search

nutofem / nuto

NuTo - yet another finite element library
https://nuto.readthedocs.io
Boost Software License 1.0
17 stars 5 forks source link

Implicit dependency of Interpolation and Intergration #34

Open TTitscher opened 7 years ago

TTitscher commented 7 years ago

eShapeType defines the shape of an InterpolationType (TRUSS, BRICK, ...) and the eTypeOrder defines the actual interpolation type and its order (EQUIDISTANT1, LOBATTO3). Independently, you have to pick integration types basically containing the same information hidden in another enum eIntegrationType, e.g. eIntegrationType::IntegrationType2D4NGauss4Ip

Problem

Solution (idea)

Maybe like:

auto triangle = NuTo::Shape::Triangle::Create()
element.SetInterpolationType(triangle.Interpolation<Equidistant>(2));
element.SetIntergrationType(triangle.Integration<Equidistant>(2));
//or
element.SetBoth<Equidistant>(triangle, 2);

Urgency

Meh. If we find a concept, this could be implemented during the PDE change... :)

Discussion

Do you see these problems as well? How should the API look like? How should the classes look like?

// template everything
NuTo::Interpolation<Triangle, Equidistant,1> 
// template some
NuTo::InterpolationTriangle<Equidistant,1>
NuTo::InterpolationTriangleEquidistant<1> 
// some dynamic
NuTo::InterpolationTriangle(EQUIDISTANT, 1)
// ...

What features would you like to see?

joergfunger commented 7 years ago

In my opinion, there is a clear distinction between integration type and interpolation type, in particular when moving to the new PDE concept. In this case, it strongly depends on the definition of the PDE (independent of the interpolation type), thus this should remain decoupled. A possible interface (instead of defining the Lobatto/Gauss) integration with the number of points would be to define the type (Gauss/Lobatto) and then the integration order, the corresponding integration type can then be chosen automatically, e.g. by also defining an order of the PDE relative to the interpolation type.

Jörg F.

2017-06-02 9:28 GMT+02:00 TTitscher notifications@github.com:

eShapeType defines the shape of an InterpolationType (TRUSS, BRICK, ...) and the eTypeOrder defines the actual interpolation type and its order ( EQUIDISTANT1, LOBATTO3). Independently, you have to pick integration types basically containing the same information hidden in another enum eIntegrationType, e.g. eIntegrationType::IntegrationType2D4NGauss4Ip

Problem

  • redundant information
  • Combination of (e.g.) LOBATTO and TETRAHEDRON can be defined, but is impossible
    • many unnecessary run time checks
  • Which integration type can accurately integrate triangles of order 4? 6IP, 10IP, 12IP?, 15IP?
    • error prone
  • matching interpolation types and integration types have to be defined in the same natural coordinate system

Solution (idea)

-

Shape classes that defines the natural coordinate system and its bounds

  • Quad - xi, eta in [-1, 1]
    • Triangle - xi, eta in triange [0,0] -- [1,0] -- [0,1];
    • ....

  • Define the general interpolation / integration concept with an enum containing e.g. Equidistant, Lobatto, Serendipity

  • Somehow create the integration types and interpolation types from an integer value containing the order

Maybe like:

auto triangle = NuTo::Shape::Triangle::Create() element.SetInterpolationType(triangle.Interpolation(2)); element.SetIntergrationType(triangle.Integration(2));//or element.SetBoth(triangle, 2);

Urgency

Meh. If we find a concept, this could be implemented during the PDE change... :) Discussion

Do you see these problems as well? How should the API look like? How should the classes look like?

// template everything NuTo::Interpolation<Triangle, Equidistant,1> // template some NuTo::InterpolationTriangle<Equidistant,1> NuTo::InterpolationTriangleEquidistant<1> // some dynamicNuTo::InterpolationTriangle(EQUIDISTANT, 1)// ...

What features would you like to see?

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/nutofem/nuto/issues/34, or mute the thread https://github.com/notifications/unsubscribe-auth/ATjSJTppIBvV5hJ5H9mKsqBRbYmhEe7tks5r_7mggaJpZM4Nt7OZ .

TTitscher commented 7 years ago

@joergfunger Could you briefly elaborate on why you prefer the decoupling? I am afraid that I overlook some issues with that.

Maybe to clarify, I was just asking myself:

To my knowledge, no (... I might be wrong in the second question, though??). And it would be nice, if these invalid choices are somehow forbidden. At compile time.

Defining the default integration order of a PDE is a neat idea IMO.

joergfunger commented 7 years ago

The check at compile time to have identical shapes (triangle, quad,..) is a very good idea, similar with the definition of an identical coordinate system. In my opinion, the only problem might be related to the order. There are several reasons, why this should be separated

2017-06-06 16:05 GMT+02:00 TTitscher notifications@github.com:

@joergfunger https://github.com/joergfunger Could you briefly elaborate on why you prefer the decoupling? I am afraid that I overlook some issues with that.

Maybe to clarify, I was just asking myself:

  • Is it meaningful to integrate a 2DTriangle interpolation with a 3DBrick integration type?
  • Is it meaningful to integrate a Lobatto interpolation with Gauss integration points?
  • Is it meaningful to perform the interpolation and the integration in different natural coordinate systems?

To my knowledge, no (... I might be wrong in the second question, though??). And it would be nice, if these invalid choices are somehow forbidden. At compile time.

Defining the default integration order of a PDE is a neat idea IMO.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/nutofem/nuto/issues/34#issuecomment-306496664, or mute the thread https://github.com/notifications/unsubscribe-auth/ATjSJbwyPPOZ6lq-zd6bO9hd0eTw8lJ7ks5sBVyqgaJpZM4Nt7OZ .

Psirus commented 6 years ago

So, here is the basic idea of Group Shape (as far as I understood it):

The shapes will be very small (ideally "member-free") classes that represent geometric shapes, allowing to match the shape of integration and interpolation for example. For now, they'll probably just have a couple of functions, for example:

class Shape
{
public:
    virtual int GeometricDimension() const = 0;
    virtual bool IsWithinShape(Eigen::VectorXd xi) const = 0; 
   //            ^ checks e.g. if xi for line is in [-1, 1]
};

We're proposing to save these by value at the relevant classes, i.e. an InterpolationBrickLinear would get a Hexahedron member and GetShape() function.

Now, the tricky bit: We thought it essential to have the "unrelated" functionality of

not included in the shape class. If you want create an interpolation from a shape, say

Interpolation& CreateInterpolation(Shape& shape, int order);

you don't have the derived class information anymore (it's just a shape), so you either need to resort to RTTI (frowned upon, for good reason), or use the Visitor pattern. With this technique (using double dispatch), we can recover the type of the shape and act accordingly. So what does it look like?

The shape gets extended to accept a Visitor class (the visitors will be our integration and interpolation type factories, etc.)

class Shape
{
public:
    virtual void accept(Visitor& visitor) const;
    // ...
};

and the visitor gets a visit function for each derived shape

class Visitor
{
public:
    virtual void visit(const Triangle&) = 0;
    virtual void visit(const Hexahedron&) = 0;
    // ...
};

Now, any derived visitor just needs to call accept on the given shape, and will "see" which of its visit functions was called.

And because I know this is very confusing in prose, take a look at a prototype interpolation factory I've implemented in the shape branch.

Any objections? Suggestions?

TTitscher commented 6 years ago

I like the idea of the shape class that may have some more geometrical information. And I like the idea that the actual creation of the interpolations/integrations is done in factories.

However, on first glance I thought of a TriangleFactory and a QuadFactory instead of a InterpolationFactory and VisualizeFactory. Why? Everything that belongs to a triangle is collected in the triangle factory. Other words: if I want to interpolate and integrate triangles, I only include the interpolation and integration for triangles. Fair enough. With the Interpolation/IntegrationFactory idea, I include every interpolation and every integration. Meh.

Another advantage is the extensibility. I am allowed to define a custom shape in my main file. Well, arguable. But this shape may come with a custom factories. This can be useful for playing around with experimental interpolations/integrations while utilizing the comfortable interface the shapes provide.

An argument I heard is that the shapes should not be extended. May be fine as well. I this case, how is the visitor pattern with its double dispatch better as a enum eShape. Then, Interpolation::GetShape() returns this enum and the InterpolationFactory accepts it. Edit: The additional features could then live in a other quasi factories like int GeometricDimension(eShape) or bool IsWithin(eShape, coords).

Psirus commented 6 years ago

TriangleFactory vs InterpolationFactory

So there are two "orthogonal" things here, the shapes themselves, and what to do with them.

First the shapes: the common shapes that come up in FEM are line, triangle, quad, tet and hex (and maybe point). With these, you have 95% of all use cases covered. VTK and Gmsh additionally support prisms and pyramids, raising that to probably 99.9%. If you still need another shape, I think it would not be unreasonable in this case to add the shape in the library itself. (Since you're already writing your own meshing and postprocessing :))

Things one can do with a shape: currently this would include creating interpolations, integration types and visualization geometries. But this list is clearly not as well defined, and what might be needed in the future here I really can't predict. Therefore, I would place more emphasis on the extensibility of the "things to be done with a shape".

Yes, the interpolation factory would need to include all interpolations that I want to be able to create. Such is life when writing factories. Conversely, a triangle factory would need to include InterpolationSimple, IntegrationTypeBase, VoronoiGeometry, AverageGeometry. And if a new "thing-that-uses-shapes" is needed, you'll have to touch all the shapes. This would the worse case in terms of difficulty to extend imho.

How is it better than an enum?

The additional features could then live in a other quasi factories like int GeometricDimension(eShape) or bool IsWithin(eShape, coords).

Alternatively, we could group these functionality by shape, and maybe collect them in some kind of entity, and maybe call these entities classes ;)

I don't think it needs to stop at these two functions, I could imagine it being useful to ask the shape "What is the shape of your faces? How many vertices do you have? What is the normal direction of your third face? Etc.". This sort of information could be useful for example in Neumann BCs or DG approaches, where one needs to integrate over the face of an element. Then, an enum + free functions becomes cumbersome, I would think.

TTitscher commented 6 years ago

My argument for the TriangleFactory was more like that it enables adding a whole new factory from the outside. Like adding TriangleDiscontinuous or whatever. But this may belong directly into NuToCore. So I think I can live with the InterpolationFactory and its siblings. In any case we should make sure that no main method requires those. For example: 

// we can now to this, which is super nice:
Cell(element)
   : mElement(element)
   , mIntegrationType(Factory(element.Interpolation())
  ...

// but we _must_ keep this to allow extensions
Cell(element, integration)
   : mElement(element)
   , mIntegration(integration)
  ...

Sorry, if this is totally obvious.

enum vs visitor: I admit that bool IsWithin(eShape, coords) is crap. But I think we will end up with easier code, if we e.g. replace Shape::accept(...) with enum Shape::Enum(). And the factories have switches like switch(shape.Enum() case TRIANGLE: .... This is like implementing a typeid by hand. So we may rely on that.

Still, we will end up with a switch in each factory that bases its behavior on the type of the shape, regardless of enum, typeinfo or visitor (see Factory::visit(Triangle&) and Factory::visit(Hex&)). Common practice would be to move that pieces of code into these classes. Just saying :cry: But just to be clear:

I think I can live with the InterpolationFactory and its siblings.

pmueller2 commented 6 years ago

I think the concept of a simple reference element of different shapes should lead to a simple object that describes just this. It has nothing to do with integration, interpolation, whatever. I usually think of interpolation on a triangle and not of a triangles interpolation facilities and this should be reflected by the code. If you use typeid or enums this will not work when using a DerivedSpecialTriangle (is probably not relevant?). The visitor pattern has no problem with that.

Psirus commented 6 years ago

So, I have continued a bit with that. The factory now uses the enum instead of the visitor, and while I have not started on the Lobatto stuff yet, I thought it would be enough to have two functions

std::unique_ptr<InterpolationSimple> CreateSerendipityInterpolation(const Shape& shape, int order);
std::unique_ptr<InterpolationSimple> CreateLobattoInterpolation(const Shape& shape, int order);

instead of introducing another enum. Would that be problematic somehow? The rest can be seen on the shapes branch.