I am importing a 3D rectangular beam mesh from Gmsh (script provided after MWE). The rectangular 3D beam has length 10, width 0.5, and height 0.5. I am defining the volume, surface, and line measures as given in the MWE below. When I assemble, I am getting the right volume. But the values of the areas of the surfaces and the lines are way off!
MWE:
from dolfin import *
import meshio
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from ufl import cofac, sqrt
from scipy import linalg as la
from ufl import *
from ufl.classes import *
import pandas as pd
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True}
import numpy
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
return out_mesh
######################################################################################################################
msh = meshio.read("mesh.msh")
######################################################################################################################
tetra_mesh = create_mesh(msh, "tetra", True)
triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", tetra_mesh)
meshio.write("surface.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh)
from dolfin import *
mesh = Mesh()
xdmf = XDMFFile(mesh.mpi_comm(),"mesh.xdmf")
xdmf.read(mesh)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
xdmf.close()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("surface.xdmf") as infile:
infile.read(mvc, "name_to_read")
sf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
xdmf.close()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-2)
with XDMFFile("mf.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds_custom = Measure("dS", domain=mesh, subdomain_data=mf) # Line measure
dx_custom = Measure("dx", domain=mesh, subdomain_data=sf) # Surface measure
dv_custom = Measure("dx", domain=mesh, subdomain_data=cf) # Volume measure
print(f"Volume = {assemble(Constant(1.0)*dv_custom(19))}") # Should be 2.5
print(f"frontbottomline = {assemble(Constant(1.0)*ds_custom(4))}") # Should be 10
print(f"fronttopline = {assemble(Constant(1.0)*ds_custom(3))}") # Should be 10
print(f"frontleftline = {assemble(Constant(1.0)*ds_custom(1))}") # Should be 0.5
print(f"frontrightline = {assemble(Constant(1.0)*ds_custom(2))}") # Should be 0.5
print(f"lefttopline = {assemble(Constant(1.0)*ds_custom(5))}") # Should be 0.5
print(f"leftbottomline = {assemble(Constant(1.0)*ds_custom(6))}") # Should be 0.5
print(f"righttopline = {assemble(Constant(1.0)*ds_custom(7))}") # Should be 0.5
print(f"rightbottomline = {assemble(Constant(1.0)*ds_custom(8))}") # Should be 0.5
print(f"backtopline = {assemble(Constant(1.0)*ds_custom(9))}") # Should be 10
print(f"backbottomline = {assemble(Constant(1.0)*ds_custom(10))}") # Should be 10
print(f"backleftline = {assemble(Constant(1.0)*ds_custom(11))}") # Should be 0.5
print(f"backleftline = {assemble(Constant(1.0)*ds_custom(12))}") # Should be 0.5
print(f"frontsurface = {assemble(Constant(1.0)*dx_custom(13))}") # Should be 5
print(f"leftsurface = {assemble(Constant(1.0)*dx_custom(14))}") # Should be 0.25
print(f"rightsurface = {assemble(Constant(1.0)*dx_custom(15))}") # Should be 0.25
print(f"topsurface = {assemble(Constant(1.0)*dx_custom(16))}") # Should be 5
print(f"bottomsurface = {assemble(Constant(1.0)*dx_custom(17))}") # Should be 5
print(f"backsurface = {assemble(Constant(1.0)*dx_custom(18))}") # Should be 5
This gives the following output:
Volume = 2.5000000000002993
frontbottomline = 0.5088177178912574
fronttopline = 0.026856850843952427
frontleftline = 0.0012651138839902972
frontrightline = 0.0012576415390319338
lefttopline = 0.23572496165217366
leftbottomline = 0.04524590534914055
righttopline = 0.003959317943007174
rightbottomline = 0.001385418482056354
backtopline = 0.02318226274902181
backbottomline = 0.023506977188747825
backleftline = 0.0009567354242365273
backleftline = 0.0012058685118405636
frontsurface = 0.04434651358155991
leftsurface = 0.002269191407108078
rightsurface = 0.0022366414437487514
topsurface = 0.04266836365558237
bottomsurface = 0.04323308955707949
backsurface = 0.04398366574788088
The values of the surface areas and length of the lines are way off!
Here is the Gmsh script:
SetFactory("OpenCASCADE");
//+
Point(1) = {0.0, 0.0, 0.0, 1.0};
//+
Point(2) = {0.0, 0.0, 0.5, 1.0};
//+
Point(3) = {10.0, 0.0, 0.5, 1.0};
//+
Point(4) = {10.0, 0.0, 0.0, 1.0};
//+
Point(5) = {0.0, 0.5, 0.0, 1.0};
//+
Point(6) = {10.0, 0.5, 0.0, 1.0};
//+
Point(7) = {10.0, 0.5, 0.5, 1.0};
//+
Point(8) = {0.0, 0.5, 0.5, 1.0};
//+
Line(1) = {8, 7};
//+
Line(2) = {7, 3};
//+
Line(3) = {3, 2};
//+
Line(4) = {2, 8};
//+
Line(5) = {5, 8};
//+
Line(6) = {1, 2};
//+
Line(7) = {1, 5};
//+
Line(8) = {5, 6};
//+
Line(9) = {6, 7};
//+
Line(10) = {1, 4};
//+
Line(11) = {4, 3};
//+
Line(12) = {4, 6};
//+
Curve Loop(1) = {1, 2, 3, 4};
//+
Plane Surface(1) = {1};
//+
Curve Loop(2) = {6, 4, -5, -7};
//+
Plane Surface(2) = {2};
//+
Curve Loop(3) = {11, -2, -9, -12};
//+
Plane Surface(3) = {3};
//+
Curve Loop(4) = {8, -12, -10, 7};
//+
Plane Surface(4) = {4};
//+
Curve Loop(5) = {11, 3, -6, 10};
//+
Plane Surface(5) = {5};
//+
Curve Loop(6) = {1, -9, -8, 5};
//+
Plane Surface(6) = {6};
//+
Surface Loop(1) = {6, 1, 3, 5, 2, 4};
//+
Volume(1) = {1};
//+
Physical Curve("fronttop", 3) = {1};
//+
Physical Curve("frontbottom", 4) = {3};
//+
Physical Curve("frontleft", 1) = {4};
//+
Physical Curve("frontright", 2) = {2};
//+
Physical Curve("lefttop", 5) = {5};
//+
Physical Curve("leftbottom", 6) = {6};
//+
Physical Curve("righttop", 7) = {9};
//+
Physical Curve("rightbottom", 8) = {11};
//+
Physical Curve("backtop", 9) = {8};
//+
Physical Curve("backbottom", 10) = {10};
//+
Physical Curve("backleft", 11) = {7};
//+
Physical Curve("backright", 12) = {12};
//+
Physical Surface("front", 13) = {1};
//+
Physical Surface("left", 14) = {2};
//+
Physical Surface("right", 15) = {3};
//+
Physical Surface("top", 16) = {6};
//+
Physical Surface("bottom", 17) = {5};
//+
Physical Surface("back", 18) = {4};
//+
Physical Volume("vol", 19) = {1};
I am not sure what I am doing wrong in defining the measures. Any suggestions would be greatly appreciated! Thanks in advance!