6
$\begingroup$

In this simplified example I consider a small FEM-mesh

<< NDSolve`FEM`
poly = Polygon[{{0, 0}, {3/2, 0}, {3/2, 1}, {1 , 1}, {1, 1/2}, {0,1/2}}]
mesh = ToElementMesh[poly, MaxCellMeasure -> 1/5,"MeshElementType" -> "TriangleElement","MeshOrder" -> 1]
mesh["Wireframe"]

enter image description here

For this mesh, I define, using ElementMeshInterpolation[], some basisfunctions φ

p = mesh ["Coordinates"] 
φ = Map[ElementMeshInterpolation[mesh, #] &,IdentityMatrix[Length[p ]]];

Trying to integrate products (matrix) of these functions

φφ =Outer[Times, Map[#[x, y] &, φ ], Map[#[x, y] &, φ ]];
M = NIntegrate[φφ, Element[{x, y}, mesh]]; // AbsoluteTiming
(*~6 seconds *)

I notice a large evaluation time.

My question: How is it possible to speed up integration of expressions formed by ElementMeshInterpolation[]-functions?

Thanks!

$\endgroup$

1 Answer 1

5
$\begingroup$
$Version

"13.0.0 for Mac OS X ARM (64-bit) (December 3, 2021)"

Many thanks to @Ulrich Neumann for the comments.

Here we provide three alternatives:

  1. Options in NIntegrate

  2. Changing the InterpolationOrder

  3. Chaning the MaxCellMeasure

In each case we provide timings and comparisons of the numerical integration.


Final version:


  • Performing the code of the OP for comparison

We have

<< NDSolve`FEM`
poly = Polygon[{{0, 0}, {3/2, 0}, {3/2, 1}, {1, 1}, {1, 1/2}, {0, 
     1/2}}];
mesh = ToElementMesh[poly, MaxCellMeasure -> 1/5, 
   "MeshElementType" -> "TriangleElement", "MeshOrder" -> 1];
mesh["Wireframe"]

mesh1

p = mesh["Coordinates"];
φ = 
  Map[ElementMeshInterpolation[mesh, #] &, 
   IdentityMatrix[Length[p]]];
φφ = 
  Outer[Times, Map[#[x, y] &, φ], 
   Map[#[x, y] &, φ]];

Ιntegrating

MOP = NIntegrate[φφ, 
    Element[{x, y}, mesh]]; // AbsoluteTiming

{2.14101, Null}


  • A suggestion at the level of NIntegrate

With the mesh of the OP untouched

<< NDSolve`FEM`
poly = Polygon[{{0, 0}, {3/2, 0}, {3/2, 1}, {1, 1}, {1, 1/2}, {0, 
     1/2}}];
mesh = ToElementMesh[poly, MaxCellMeasure -> 1/5, 
   "MeshElementType" -> "TriangleElement", "MeshOrder" -> 1];
mesh["Wireframe"]

the following

p = mesh["Coordinates"]
φ = 
  Map[ElementMeshInterpolation[mesh, #] &, 
   IdentityMatrix[Length[p]]];
φφ = 
  Outer[Times, Map[#[x, y] &, φ], 
   Map[#[x, y] &, φ]];
M1 = NIntegrate[φφ, Element[{x, y}, mesh], 
    Method -> {"FiniteElement", 
      "MeshOptions" -> {"MeshElementType" -> TetrahedronElement, 
        "MaxCellMeasure" -> 1/5}}]; // AbsoluteTiming

returns

{0.542303, Null}

  • Comparing the numerical results

A simple comparison

MOP - M1 // Chop[#, 7 10^-3] &

compare1


  • A suggestion at the level of ElementMeshInterpolation

With the mesh of the OP untouched

<< NDSolve`FEM`
poly = Polygon[{{0, 0}, {3/2, 0}, {3/2, 1}, {1, 1}, {1, 1/2}, {0, 
     1/2}}];
mesh = ToElementMesh[poly, MaxCellMeasure -> 1/5, 
   "MeshElementType" -> "TriangleElement", "MeshOrder" -> 1];
mesh["Wireframe"]

we do

p = mesh["Coordinates"]
φ = 
  Map[ElementMeshInterpolation[mesh, #, InterpolationOrder -> 1] &, 
   IdentityMatrix[Length[p]]];
φφ = 
  Outer[Times, Map[#[x, y] &, φ], 
   Map[#[x, y] &, \[CurlyPhi]]];
M2 = NIntegrate[φφ, 
    Element[{x, y}, mesh]]; // AbsoluteTiming

to get

{2.13627, Null}

which is not a huge improvement.

  • Comparing the numerical results

As before

MOP - M2 

compare2


  • A suggestion at the level of MaxCellMeasure

We begin with:

<< NDSolve`FEM`
poly = Polygon[{{0, 0}, {3/2, 0}, {3/2, 1}, {1, 1}, {1, 1/2}, {0, 
     1/2}}];
mesh = ToElementMesh[poly, MaxCellMeasure -> 1/4, 
   "MeshElementType" -> "TriangleElement", "MeshOrder" -> 1];
mesh["Wireframe"]

mesh2

and proceed as before

p = mesh["Coordinates"]
φ = 
  Map[ElementMeshInterpolation[mesh, #] &, IdentityMatrix[Length[p]]];

φφ = 
  Outer[Times, Map[#[x, y] &, φ], 
   Map[#[x, y] &, φ]];
M3 = NIntegrate[φφ, 
    Element[{x, y}, mesh]]; // AbsoluteTiming

{1.71267, Null}

  • Comparing the numerical results

As before

MOP[[All ;; -2, All ;; -2]] - M3 // Chop[#, 10^-3] &

compare3


$\endgroup$
8
  • 1
    $\begingroup$ Thank you for your answer. Evaluation of your code gives step1: 5.84s (9points) and step2: 4.72s (8points). Roughly spoken the speedup depends, not surprising, on the number of points. Still I'm looking for a speedup which is comparable to FEM realization in Mathematica $\endgroup$ Mar 24 at 7:12
  • $\begingroup$ @UlrichNeumann thanks for taking the time to check. I am a bit confused that the first integration gives roughly 6s, because I copied exactly what you wrote in the OP. And it is evident that I misunderstood your needs, sorry about that. $\endgroup$
    – bmf
    Mar 24 at 7:14
  • $\begingroup$ @UlrichNeumann sorry for the double-comment, but I forgot to say that if there are no follow-ups here, I feel I should delete my response as it is not well-suited. So just let me know if you have additional comments. Thanks :-) $\endgroup$
    – bmf
    Mar 24 at 7:17
  • 1
    $\begingroup$ Thank you for your effort. I checked your and my code and corrected my answer to 6s evaluation time. $\endgroup$ Mar 24 at 8:31
  • 1
    $\begingroup$ Thank you for your helpful suggestions. I didn't know NIntegrate-Method-FiniteElement, seems to factor 10 faster! $\endgroup$ Mar 24 at 22:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.