Hi aherrema
For your mesh question:- if you want to use old dolfin then you can use a mixed dimensional mesh with MeshView or a 2d mesh in a 3d space. I personally wouldn’t bother with quad meshes, get it working with tris first.
Here is some sample code I used to make the mesh for my project - it won’t work directly for you because it relies on other functions as well but you can get the idea from it
Note the lines which create a VectorFunctionSpace and dolfin.File(self.Name+’.xml.gz’) <<self._dolfin_mesh are important - they are doing nothing but for some reason the code would crash without them
def dolfin_meshes(self, *, meshType='mixedDimensional',gdim=3):
cellTypes={1:df.CellType.Type.interval.__repr__().split('.')[-1] , 2: df.CellType.Type.triangle.__repr__().split('.')[-1] ,3:df.CellType.Type.tetrahedron.__repr__().split('.')[-1]}
if 'separate' in meshType.lower():
ret=[None]*3
for dim in [1, 2, 3]:
elts, nodes=self.elts_and_nodes_of_dim(dim)
if len(elts):
model=self.tri_model(only2d=1) if dim==2 else self
gdim1= dim if 'dims' in meshType.lower() else gdim
mesh = df.Mesh()
mesh_editor = df.MeshEditor()
mesh_editor.open(mesh,cellTypes[dim], dim, gdim1)
mesh_editor.init_vertices(len(nodes))
mesh_editor.init_cells(len(elts ))#
nodedict={nd.Id:i for i, nd in enumerate(nodes)}
for i, nd in enumerate(nodes):mesh_editor.add_vertex(i,nd.x[:gdim1])
for i, el in enumerate(elts): mesh_editor.add_cell(i,numpy.array([nodedict[n.Id] for n in el.nodes], dtype=numpy.uintp))
mesh_editor.close()
VectorSpace1= df.VectorFunctionSpace(mesh, 'Lagrange', 1, dim=gdim)#for some reason this seems needed to make fenics initiate the facets
ret[dim-1]=mesh
return ret
if 'mixed' in meshType.lower():
if self._dolfin_mesh==None:
self._dolfin_mesh = df.Mesh()
mesh_editor = df.MeshEditor()
gdim=3
cellTypes={1:df.CellType.Type.interval.__repr__().split('.')[-1] , 2: df.CellType.Type.triangle.__repr__().split('.')[-1] ,3:df.CellType.Type.tetrahedron.__repr__().split('.')[-1]}
dim=2 if len(self.shell_elts()[0]) else 1 #TODO deal with adding different dims
model=self
maxdim=3-[len(model.elts_of_dim(dim)[0])>0 for dim in [1, 2, 3]][::-1].index(1)
if maxdim==2:
elts=self.tri_model(only2d=1).elts_of_dim(maxdim)[0]
else:
elts=model.elts_of_dim(maxdim)[0]
mesh_editor.open(self._dolfin_mesh,cellTypes[maxdim], maxdim, gdim)
mesh_editor.init_vertices(len(model.nodes))
mesh_editor.init_cells(len(elts ))#TODO this is probably wrong for mixed dim meshes
for i, nd in enumerate(model.nodes):mesh_editor.add_vertex(i,nd.x)
if len(elts):
for i, el in enumerate(elts):
mesh_editor.add_cell(i,numpy.array([model.nodedict[n.Id] for n in el.nodes], dtype=numpy.uintp))
mesh_editor.close()
dolfin.File(self.Name+'.xml.gz') <<self._dolfin_mesh
self._dolfin_meshes =[None]*3
for dim in [1, 2, 3]:
if dim==maxdim:
self._dolfin_meshes[dim-1]=self._dolfin_mesh
elif len(self.elts_of_dim(dim)[0]) :
marker =dolfin.MeshFunction("size_t", self._dolfin_mesh, dim, 0) # mesh.topology().dim() - 1, 0)
for f in dolfin.facets(self._dolfin_mesh): marker[f] = 1
self._dolfin_meshes[dim-1] = dolfin.MeshView.create(marker, 1)#dim)
return self._dolfin_meshes
Jermey Bleyer is right about fenics shells (its 3D example is a projection to 2D) but you could extend fenics shells with the same approach Jeremy used for beams, I did that a couple of years ago (i.e. use ideas from Elastic 3D beam structures — Numerical tours of continuum mechanics using FEniCS master documentation to update the fenics shells example. It worked in principle, even for large deformations, but it had some strange issues when the node order of adjacent elements changed , my maths was not strong enough to resolve it.
Using Jeremy’s new shell code will probably be a better starting point - his 3D beam code is well laid out, its easy to adapt his 3D beam example to work an arbitrary beam mesh, I guess his new shell code will also be clear. You could probably add the fenics shell implementations to it later. If you do add the fenics shells then I would recommend to start with the PSRI element, it performs well and is fairly easy to understand, the MITC elements are a bit more complex.
There is some IGA based code from David Kamensky which handles 3d shells - tIGAr , I haven’t used it myself but the code looks good and that also might be useful to you.