Hi
I have used FEniCS 2019 version to solve reaction-diffusion PDEs on immersed manifolds e.g surfaces embedded in 3D Euclidean space. The paper that described this feature of FEniCS is Automating the solution of PDEs on the sphere and other manifolds in FEniCS 1.2.
Is this feature available in dolfinx already?
Regards,
As far as Im aware these features are in DOLFINx.
For instance, you can adapt Figure 3 to work with ffcx:
from ufl import Cell, VectorElement, FunctionSpace, Mesh, Coefficient
# Define triangle cell embedded in Rˆ3
cell = Cell("triangle", 3)
# Define Lagrange vector element
coordinate_element = VectorElement("Lagrange", cell, 1)
mesh = Mesh(coordinate_element)
V = FunctionSpace(mesh, coordinate_element)
# Arguments defined over V will have 3 components:
u = Coefficient(V)
print(u[0], u[1], u[2])
python3 -m ffcx manifold.py
w_0[0] w_0[1] w_0[2]
or Figure 4:
from ufl import Cell, VectorElement, FunctionSpace, Mesh, Coefficient, FiniteElement, TrialFunction, TestFunction, inner, grad, dx
# Define triangle cell embedded in Rˆ3
cell = Cell("triangle", 3)
# Define Lagrange vector element
coordinate_element = VectorElement("Lagrange", cell, 1)
mesh = Mesh(coordinate_element)
element = FiniteElement("Lagrange", cell, 1)
V = FunctionSpace(mesh, element)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v))*dx
which yields the kernel in manifold.c
void tabulate_tensor_integral_6a56ba35a29fca58aae87b107061fd2de6b29fd1(double* restrict A,
const double* restrict w,
const double* restrict c,
const double* restrict coordinate_dofs,
const int* restrict entity_local_index,
const uint8_t* restrict quadrature_permutation)
{
// Quadrature rules
static const double weights_083[1] = { 0.5 };
// Precomputed values of basis functions and precomputations
// FE* dimensions: [permutation][entities][points][dofs]
static const double FE10_C0_D010_Q083[1][1][1][3] = { { { { -1.0, 0.0, 1.0 } } } };
static const double FE6_C0_D100_Q083[1][1][1][3] = { { { { -1.0, 1.0, 0.0 } } } };
// Quadrature loop independent computations for quadrature rule 083
double J_c1 = 0.0;
double J_c0 = 0.0;
double J_c2 = 0.0;
double J_c4 = 0.0;
double J_c3 = 0.0;
double J_c5 = 0.0;
for (int ic = 0; ic < 3; ++ic)
{
J_c1 += coordinate_dofs[ic * 3] * FE10_C0_D010_Q083[0][0][0][ic];
J_c0 += coordinate_dofs[ic * 3] * FE6_C0_D100_Q083[0][0][0][ic];
J_c2 += coordinate_dofs[ic * 3 + 1] * FE6_C0_D100_Q083[0][0][0][ic];
J_c4 += coordinate_dofs[ic * 3 + 2] * FE6_C0_D100_Q083[0][0][0][ic];
J_c3 += coordinate_dofs[ic * 3 + 1] * FE10_C0_D010_Q083[0][0][0][ic];
J_c5 += coordinate_dofs[ic * 3 + 2] * FE10_C0_D010_Q083[0][0][0][ic];
}
double sp_083[74];
sp_083[0] = J_c0 * J_c0;
sp_083[1] = J_c2 * J_c2;
sp_083[2] = sp_083[0] + sp_083[1];
sp_083[3] = J_c4 * J_c4;
sp_083[4] = sp_083[2] + sp_083[3];
sp_083[5] = J_c1 * J_c1;
sp_083[6] = J_c3 * J_c3;
sp_083[7] = sp_083[5] + sp_083[6];
sp_083[8] = J_c5 * J_c5;
sp_083[9] = sp_083[7] + sp_083[8];
sp_083[10] = sp_083[4] * sp_083[9];
sp_083[11] = J_c0 * J_c1;
sp_083[12] = J_c2 * J_c3;
sp_083[13] = sp_083[11] + sp_083[12];
sp_083[14] = J_c4 * J_c5;
sp_083[15] = sp_083[13] + sp_083[14];
sp_083[16] = sp_083[15] * sp_083[15];
sp_083[17] = sp_083[10] + -1 * sp_083[16];
sp_083[18] = sp_083[4] / sp_083[17];
sp_083[19] = J_c1 * sp_083[18];
sp_083[20] = (-1 * sp_083[15]) / sp_083[17];
sp_083[21] = J_c0 * sp_083[20];
sp_083[22] = sp_083[19] + sp_083[21];
sp_083[23] = sp_083[9] / sp_083[17];
sp_083[24] = J_c0 * sp_083[23];
sp_083[25] = J_c1 * sp_083[20];
sp_083[26] = sp_083[24] + sp_083[25];
sp_083[27] = sp_083[22] * sp_083[22];
sp_083[28] = sp_083[22] * sp_083[26];
sp_083[29] = sp_083[26] * sp_083[26];
sp_083[30] = J_c3 * sp_083[18];
sp_083[31] = J_c2 * sp_083[20];
sp_083[32] = sp_083[30] + sp_083[31];
sp_083[33] = J_c2 * sp_083[23];
sp_083[34] = J_c3 * sp_083[20];
sp_083[35] = sp_083[33] + sp_083[34];
sp_083[36] = sp_083[32] * sp_083[32];
sp_083[37] = sp_083[32] * sp_083[35];
sp_083[38] = sp_083[35] * sp_083[35];
sp_083[39] = sp_083[27] + sp_083[36];
sp_083[40] = sp_083[28] + sp_083[37];
sp_083[41] = sp_083[29] + sp_083[38];
sp_083[42] = J_c5 * sp_083[18];
sp_083[43] = J_c4 * sp_083[20];
sp_083[44] = sp_083[42] + sp_083[43];
sp_083[45] = J_c4 * sp_083[23];
sp_083[46] = J_c5 * sp_083[20];
sp_083[47] = sp_083[45] + sp_083[46];
sp_083[48] = sp_083[44] * sp_083[44];
sp_083[49] = sp_083[44] * sp_083[47];
sp_083[50] = sp_083[47] * sp_083[47];
sp_083[51] = sp_083[39] + sp_083[48];
sp_083[52] = sp_083[40] + sp_083[49];
sp_083[53] = sp_083[41] + sp_083[50];
sp_083[54] = J_c2 * J_c5;
sp_083[55] = J_c3 * J_c4;
sp_083[56] = sp_083[54] + -1 * sp_083[55];
sp_083[57] = sp_083[56] * sp_083[56];
sp_083[58] = J_c1 * J_c4;
sp_083[59] = J_c0 * J_c5;
sp_083[60] = sp_083[58] + -1 * sp_083[59];
sp_083[61] = sp_083[60] * sp_083[60];
sp_083[62] = sp_083[57] + sp_083[61];
sp_083[63] = J_c0 * J_c3;
sp_083[64] = J_c1 * J_c2;
sp_083[65] = sp_083[63] + -1 * sp_083[64];
sp_083[66] = sp_083[65] * sp_083[65];
sp_083[67] = sp_083[62] + sp_083[66];
sp_083[68] = sqrt(sp_083[67]);
sp_083[69] = 1 * sp_083[68];
sp_083[70] = fabs(sp_083[69]);
sp_083[71] = sp_083[51] * sp_083[70];
sp_083[72] = sp_083[52] * sp_083[70];
sp_083[73] = sp_083[53] * sp_083[70];
for (int iq = 0; iq < 1; ++iq)
{
const double fw0 = sp_083[73] * weights_083[iq];
const double fw1 = sp_083[72] * weights_083[iq];
const double fw2 = sp_083[71] * weights_083[iq];
double t0[3];
double t1[3];
for (int i = 0; i < 3; ++i)
{
t0[i] = fw0 * FE6_C0_D100_Q083[0][0][0][i] + fw1 * FE10_C0_D010_Q083[0][0][0][i];
t1[i] = fw1 * FE6_C0_D100_Q083[0][0][0][i] + fw2 * FE10_C0_D010_Q083[0][0][0][i];
}
for (int i = 0; i < 3; ++i)
for (int j = 0; j < 3; ++j)
A[3 * i + j] += FE6_C0_D100_Q083[0][0][0][j] * t0[i] + FE10_C0_D010_Q083[0][0][0][j] * t1[i];
}
}
which is different from using
cell = Cell("triangle", 2)
yielding
void tabulate_tensor_integral_27b24fdf56699cf003a569e3622e9eaddc6a64cf(double* restrict A,
const double* restrict w,
const double* restrict c,
const double* restrict coordinate_dofs,
const int* restrict entity_local_index,
const uint8_t* restrict quadrature_permutation)
{
// Quadrature rules
static const double weights_083[1] = { 0.5 };
// Precomputed values of basis functions and precomputations
// FE* dimensions: [permutation][entities][points][dofs]
static const double FE4_C0_D10_Q083[1][1][1][3] = { { { { -1.0, 1.0, 0.0 } } } };
static const double FE8_C0_D01_Q083[1][1][1][3] = { { { { -1.0, 0.0, 1.0 } } } };
// Quadrature loop independent computations for quadrature rule 083
double J_c0 = 0.0;
double J_c3 = 0.0;
double J_c1 = 0.0;
double J_c2 = 0.0;
for (int ic = 0; ic < 3; ++ic)
{
J_c0 += coordinate_dofs[ic * 3] * FE4_C0_D10_Q083[0][0][0][ic];
J_c3 += coordinate_dofs[ic * 3 + 1] * FE8_C0_D01_Q083[0][0][0][ic];
J_c1 += coordinate_dofs[ic * 3] * FE8_C0_D01_Q083[0][0][0][ic];
J_c2 += coordinate_dofs[ic * 3 + 1] * FE4_C0_D10_Q083[0][0][0][ic];
}
double sp_083[20];
sp_083[0] = J_c0 * J_c3;
sp_083[1] = J_c1 * J_c2;
sp_083[2] = sp_083[0] + -1 * sp_083[1];
sp_083[3] = J_c0 / sp_083[2];
sp_083[4] = (-1 * J_c1) / sp_083[2];
sp_083[5] = sp_083[3] * sp_083[3];
sp_083[6] = sp_083[3] * sp_083[4];
sp_083[7] = sp_083[4] * sp_083[4];
sp_083[8] = J_c3 / sp_083[2];
sp_083[9] = (-1 * J_c2) / sp_083[2];
sp_083[10] = sp_083[9] * sp_083[9];
sp_083[11] = sp_083[8] * sp_083[9];
sp_083[12] = sp_083[8] * sp_083[8];
sp_083[13] = sp_083[5] + sp_083[10];
sp_083[14] = sp_083[6] + sp_083[11];
sp_083[15] = sp_083[12] + sp_083[7];
sp_083[16] = fabs(sp_083[2]);
sp_083[17] = sp_083[13] * sp_083[16];
sp_083[18] = sp_083[14] * sp_083[16];
sp_083[19] = sp_083[15] * sp_083[16];
for (int iq = 0; iq < 1; ++iq)
{
const double fw0 = sp_083[19] * weights_083[iq];
const double fw1 = sp_083[18] * weights_083[iq];
const double fw2 = sp_083[17] * weights_083[iq];
double t0[3];
double t1[3];
for (int i = 0; i < 3; ++i)
{
t0[i] = fw0 * FE4_C0_D10_Q083[0][0][0][i] + fw1 * FE8_C0_D01_Q083[0][0][0][i];
t1[i] = fw1 * FE4_C0_D10_Q083[0][0][0][i] + fw2 * FE8_C0_D01_Q083[0][0][0][i];
}
for (int i = 0; i < 3; ++i)
for (int j = 0; j < 3; ++j)
A[3 * i + j] += FE4_C0_D10_Q083[0][0][0][j] * t0[i] + FE8_C0_D01_Q083[0][0][0][j] * t1[i];
}
}
Thank you very much for the prompt response, @dokken! Thank you for sharing the generated code as evidence.