Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Redo the quadrilateral measure function. #15755

Merged
merged 1 commit into from Jul 18, 2023

Conversation

drwells
Copy link
Member

@drwells drwells commented Jul 15, 2023

Supersedes #10419.

This is a lot simpler and somewhat more efficient:

  • master, best of 10, summing hyper ball face measures: 2434 ms
  • feature, same setup: 2124 ms

Lets not use complex checks for planarity - either we have a parallelogram (up to the last few significant digits) or we do not. If people are reading in CAD data then they probably don't have uniform meshes anyway. I don't think we should lower the tolerance in this case either: we should compute the most accurate answer we can.

benchmark source:

#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria.h>

#include <chrono>
#include <iostream>

int
main()
{
  using namespace dealii;

  Triangulation<3> tria;
  GridGenerator::hyper_ball(tria);

  tria.refine_global(3);

  double area = 0.0;
  auto   t0   = std::chrono::high_resolution_clock::now();
  for (unsigned int i = 0; i < 1000; ++i)
    for (const auto &face : tria.active_face_iterators())
      area += face->measure();
  auto t1 = std::chrono::high_resolution_clock::now();
  std::cout
    << "time elapsed = "
    << std::chrono::duration_cast<std::chrono::milliseconds>(t1 - t0).count()
    << std::endl;
}

Copy link
Member

@marcfehling marcfehling left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks way simpler, I like it!

source/grid/tria_accessor.cc Outdated Show resolved Hide resolved
source/grid/tria_accessor.cc Outdated Show resolved Hide resolved
source/grid/tria_accessor.cc Show resolved Hide resolved
source/grid/tria_accessor.cc Show resolved Hide resolved
@kronbichler
Copy link
Member

kronbichler commented Jul 17, 2023

@drwells I am curious about the order of the quadrature formula. As far as I see, we have the cross product, which contains the square of some bilinear terms, but then we take the norm and thus a square root. So the integrand is non-polynomial, but not near the singular spot for sqrt because we have positive determinants for the cells. This means that using 4 points should give a really good approximation. Is your reasoning along those points?

Copy link
Member

@kronbichler kronbichler left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I like the simpler formula. Thank you.

@drwells
Copy link
Member Author

drwells commented Jul 17, 2023

Is your reasoning along those points?

Yes - to expand a little, the normal vector is constant when a = 0 and so for faces that are not extremely distorted (if I recall correctly, Strang and Fix prove that things won't work if cells have curvature on the same order of their diameter, so in that case we have bigger problems) we are integrating a function which varies very little. Similarly, the two tangential components correspond to the two reference cell independent variables: i.e., the tangent vectors should be close to 90 degrees apart unless we have a pinched quadrilateral (and in that case we also have bad problems outside of this function) which avoids a near-zero cross product.

More simply, though, the previous version of this function used the same 16-point quadrature formula - it just unrolled all the loops manually. This version is compatible and works with the test suite. I don't want to break people's stuff by picking a different quadrature rule.

This is a lot simpler and somewhat more efficient:
master, best of 10, summing hyper ball face measures: 2434 ms
feature, same setup: 2124 ms

Lets not use complex checks for planarity - either we have a parallelogram (up
to the last few significant digits) or we do not. If people are reading in CAD
data then they probably don't have uniform meshes anyway.

P.S.

benchmark source:

#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria.h>

#include <chrono>
#include <iostream>

int
main()
{
  using namespace dealii;

  Triangulation<3> tria;
  GridGenerator::hyper_ball(tria);

  tria.refine_global(3);

  double area = 0.0;
  auto   t0   = std::chrono::high_resolution_clock::now();
  for (unsigned int i = 0; i < 1000; ++i)
    for (const auto &face : tria.active_face_iterators())
      area += face->measure();
  auto t1 = std::chrono::high_resolution_clock::now();
  std::cout
    << "time elapsed = "
    << std::chrono::duration_cast<std::chrono::milliseconds>(t1 - t0).count()
    << std::endl;
}
@tamiko tamiko merged commit 0265947 into dealii:master Jul 18, 2023
15 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

4 participants