search

mfem / mfem

Lightweight, general, scalable C++ library for finite element methods
http://mfem.org
BSD 3-Clause "New" or "Revised" License
1.74k stars 502 forks source link

Quadrature order for BoundaryNormalLFIntegrator etc. #2002

Open jandrej opened 3 years ago

jandrej commented 3 years ago

I ran into an issue where the culprit was the quadrature order of BoundaryNormalLFIntegrator. I was expecting a code to reach a certain order of convergence (in time) but at a certain point this integrator accumulated too much error and went unstable. Increasing the quadrature order fixed that.

Now the suggestion is to check all linear integrators for proper integration rules. I know those might be highly application dependent, but it might be better to have safe defaults and leave the option to reduce the order in specific cases.

tzanio commented 3 years ago

I agree -- we should discuss a consistent rule for assigning the integration order and update all integrators accordingly.

Lets have this discussion in this issue👇

There are several factors that come into play

  1. the polynomial order of the finite element space(s), fem_order
  2. the polynomial order of the geometry mapping (Jacobian), geom_order

Coefficients, can generally be anything, but using projection or just thinking of them as being approximated we can also consider

  1. the polynomial order of approximation for the coefficient, coeff_order

Based on 1-3, we can in principle determine an order for the integrand as a polynomial or a rational function. We can then proceed with

  1. defining an approximation polynomial order for the integrand, order
  2. defining a quadrature order, depending on the accuracy we want, quad_order

@mfem/developers: Is this an acceptable way to describe the factors determining the quadrature? Anything we are missing?

To start the discussion with a concrete example -- for the case of DomainLFIntegrator what do you think should be coeff_order, order and quad_order for given fem_order and geom_order?

(Currently it is quad_order = 2*fem_order, see the comments in lininteg.hpp)

stale[bot] commented 3 years ago

:warning: This issue or PR has been automatically marked as stale because it has not had any activity in the last month. If no activity occurs in the next week, it will be automatically closed. Thank you for your contributions.

marc-85 commented 3 years ago

I have an issue which might be related to this. In ex18.hpp, the integration rule for face integration is defined with const IntegrationRule *ir = &IntRules.Get(Tr.GetGeometryType(), intorder); with intorder defined earlier and intrder>0.

Now, when I try ir->GetOrder() I get 0. Am I doing something wrong? Many thanks, marc

jandrej commented 3 years ago

Are you talking about a 1D example or the geometry is a "point" by any chance? If you take a look at https://github.com/mfem/mfem/blob/master/fem/intrules.cpp#L921, the code sets the order to zero if the element type is a point.

marc-85 commented 3 years ago

I'm sorry, I didn't mention. I am running the Euler vortex (2D) with the periodic-square.mesh