Fast integration of ElementMeshInterpolation[]-functions

Question

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"]

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!

Answer 1

$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"]

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] &

---

• 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 

---

• 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"]

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] &

---