2

Context

I would like to (partially) answer my own question here (ok its a bit cheesy but...)

Question

I am interested in defining an indicator function which value would be 1 on a cell and zero outside. I am hoping to use this with the FEM package.

Example

For instance, let me define a set of 4 cells:

Needs["NDSolve`FEM`"];
reg0 = Rectangle[{0, 0}, {1, 1}];
mesh0 =  ToElementMesh[reg0, MaxCellMeasure -> 0.5, AccuracyGoal -> 0]
mesh0["Wireframe"]

Mathematica graphics

I can plot a function which changes value on each cell:

idx = mesh0["MeshElements"][[1, 1]]; 
Table[m1 = 
   ToElementMesh[mesh0["Coordinates"][[ idx[[i]]]], 
    MaxCellMeasure -> 1, AccuracyGoal -> 0]; 
  Plot3D[i, {x, y} \[Element] m1], {i, 1, Length[idx]}] // Show

Mathematica graphics

so I am not far, but What I want is to be able to build

$$ F(x,y) = 1 \quad \mbox{if} \quad {x,y} \in Cell_i $$

I am fairly certain there must be a simple elegant solution to this small problem

Constraint

I would like a solution which does not assume that the cells a necessarily squares: e.g. it should also work for

reg0 = Disk[]
mesh0 = ToElementMesh[reg0, MaxCellMeasure -> 0.5, AccuracyGoal -> 0]
mesh0["Wireframe"]

Mathematica graphics

Ideally it should also work in 3D as well.

Possible generalisation

It would be of interest to be able to define BSpline basis over such mesh element?

chris
  • 22,860
  • 5
  • 60
  • 149

1 Answers1

1

Here is my feeble attempt: there might be much more efficient methods around?

Needs["NDSolve`FEM`"];
reg0 = Disk[]
mesh0 = ToElementMesh[reg0, MaxCellMeasure -> 0.5, AccuracyGoal -> 0]
idx = mesh0["MeshElements"][[1, 1]];
pol = Table[
  mesh0["Coordinates"][[ idx[[i]]]] // ConvexHullMesh, {i, 
   Length[idx]}]

Mathematica graphics

Now I can define the indicator of the second cell as

F[x_, y_] := Boole[ RegionMember[pol[[2]], {x, y}]]

so that

Plot3D[F[x, y], {x, y} \[Element] mesh0, PlotPoints -> 30, 
 PlotTheme -> "Scientific"]

Mathematica graphics

Note that the strategy works in 3D as well

reg0 = Tetrahedron[];
mesh0 = ToElementMesh[reg0, MaxCellMeasure -> 0.5, AccuracyGoal -> 0]
pol = Table[mesh0["Coordinates"][[ idx[[i]]]] // ConvexHullMesh, {i,Length[idx]}]
idx = mesh0["MeshElements"][[1, 1]]

Mathematica graphics

F[x_, y_, z_] := Boole[ RegionMember[pol[[1]], {x, y, z}]]

Mathematica graphics

chris
  • 22,860
  • 5
  • 60
  • 149