UFL Measures on Sub-domains not Evaluated

I have a mesh imported from gmsh with defined facet tags for Surfaces and Edges. Now If I want to check area and length of this sub-domains, I am getting zero value of length of subdomains. However, for area, I am able to evaluate. Also, for the complete domain, I am able to calculate the both length and area. What am I doing wrong here.

import ufl

dx = ufl.Measure("dx", domain=domain, subdomain_data=cell_markers)
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_markers)

area = fem.assemble_scalar(fem.form(1.0 * dx(1)))
print(f"Computed area = {area:.3f}")


length_1 = fem.assemble_scalar(fem.form(1.0 * ds(1)))
print(f"Computed length = {length_1:.3f}")

To calculate the «length» of a domain, you would need an interface marker of the facets of the subdomain, or you can create a submesh for each subdomain and then evaluate ds as done for the full mesh.

Note that your example is not reproducible as you haven’t provided a mesh.

If you want to make an interface marker after meshing see

in Discontinuity at interface using mixed domains - #2 by dokken

I have defined the length markers in the mesh as follows.

cl__1 = 0.001;
Point(1) = {0.01, 0.07000000000000001, 0, cl__1};
Point(2) = {0.03, 0.07000000000000001, 0, cl__1};
Point(3) = {0.03, 0.08, 0, cl__1};
Point(4) = {0.01, 0.08, 0, cl__1};
Point(5) = {0.04, 0.07000000000000001, 0, cl__1};
Point(6) = {0.06, 0.07000000000000001, 0, cl__1};
Point(7) = {0.06, 0.08, 0, cl__1};
Point(8) = {0.04, 0.08, 0, cl__1};
Point(9) = {0.07000000000000001, 0.07000000000000001, 0, cl__1};
Point(10) = {0.09, 0.07000000000000001, 0, cl__1};
Point(11) = {0.09, 0.08, 0, cl__1};
Point(12) = {0.07000000000000001, 0.08, 0, cl__1};
Point(13) = {0, 0.06, 0, cl__1};
Point(14) = {0.1, 0.06, 0, cl__1};
Point(15) = {0.1, 0.09, 0, cl__1};
Point(16) = {0, 0.09, 0, cl__1};
Point(17) = {0, 0, 0, cl__1};
Point(18) = {0.1, 0, 0, cl__1};
Point(19) = {0.15, 0.1, 0, cl__1};
Point(20) = {-0.05, 0.1, 0, cl__1};
Line(1) = {4, 1};
Line(2) = {1, 2};
Line(3) = {2, 3};
Line(4) = {3, 4};
Line(5) = {8, 5};
Line(6) = {5, 6};
Line(7) = {6, 7};
Line(8) = {7, 8};
Line(9) = {12, 9};
Line(10) = {9, 10};
Line(11) = {10, 11};
Line(12) = {11, 12};
Line(13) = {13, 14};
Line(14) = {14, 15};
Line(15) = {15, 16};
Line(16) = {16, 13};
Line(17) = {17, 18};
Line(18) = {18, 19};
Line(19) = {19, 20};
Line(20) = {20, 17};
Curve Loop(1) = {1, 2, 3, 4};
Plane Surface(1) = {1};
Curve Loop(2) = {5, 6, 7, 8};
Plane Surface(2) = {2};
Curve Loop(3) = {9, 10, 11, 12};
Plane Surface(3) = {3};
Curve Loop(4) = {1, 2, 3, 4, -8, -7, -6, -5, -12, -11, -10, -9, -15, -14, -13, -16};
Plane Surface(4) = {4};
Curve Loop(5) = {16, 13, 14, 15, -20, -19, -18, -17};
Plane Surface(5) = {5};
Physical Line("bottom") = {17};
Physical Line("right") = {18};
Physical Line("top") = {19};
Physical Line("left") = {20};
Physical Line("EF") = {13};
Physical Line("FG") = {14};
Physical Line("GH") = {19};
Physical Line("HE") = {16};
Physical Line("IJ") = {2};
Physical Line("JK") = {3};
Physical Line("KL") = {4};
Physical Line("LI") = {1};
Physical Line("MN") = {6};
Physical Line("NO") = {7};
Physical Line("OP") = {8};
Physical Line("PM") = {5};
Physical Line("QR") = {10};
Physical Line("RS") = {11};
Physical Line("ST") = {12};
Physical Line("TQ") = {9};
Physical Surface("my_surface") = {5};
Physical Surface("EFGH") = {4};
Physical Surface("IJKL") = {1};
Physical Surface("MNOP") = {2};
Physical Surface("QRST") = {3};
Field[1] = Box;
Field[1].Thickness = 0;
Field[1].VIn = 0;
Field[1].VOut = 0;
Field[1].XMax = 0;
Field[1].XMin = 0;
Field[1].YMax = 0;
Field[1].YMin = 0;
Field[1].ZMax = 0;
Field[1].ZMin = 0;

Change this to dS then, as ds means exterior facet, ie any facet that is connected to only one cell.

1 Like