Consistent non-normalized normal field for immersed mesh?

1

Hi, I am trying to construct a consistently oriented non-normalized normal field for an immersed mesh. Below is my best attempt, trying two slightly different things. Neither give the result I need.

from dolfin import *
from ufl import Jacobian
# create immersed mesh with non-uniform Jacobian
mesh = BoxMesh(Point(0, 0, 0), Point(2, 2, 2), 2, 2, 2)
mesh.coordinates()[:,0] *= .5   
mesh.coordinates()[:,1] *= .5
srf = BoundaryMesh(mesh, "exterior", order=True)  #False --> can't project
J = Jacobian(srf)

# Tangent basis
G1 = J[:,0]
G2 = J[:,1]
N = cross(G1, G2) # expect this to be outward facing scaled normal field

# Visualize 
Nh = project(N, VectorFunctionSpace(srf, 'DG', degree=2, dim=3))
Nh.rename('n','n')
f = XDMFFile('Nh.xdmf')  # plotting in red
f.write_checkpoint(Nh,'n', 0)
f.close()

# This also didn't work (using CellNormal for orientation and norm of N for magnitude)
global_normal = Expression(("x[0]", "x[1]", "x[2]"), degree=1)
srf.init_cell_orientations(global_normal)
n = CellNormal(srf)
N1 = sqrt(dot(N,N))*n
Nh1 = project(N1, VectorFunctionSpace(srf, 'DG', degree=2, dim=3))
Nh1.rename('n','n')
f = XDMFFile('Nh1.xdmf')  # plotting in blue
f.write_checkpoint(Nh1,'n', 0)
f.close()

As you can can see in the screenshot, the second method gives slightly better results, but many of the normals still point inwards. What am I doing wrong? Thanks in advance!

image
image772Ă—1055 164 KB

2

Hi @niewiarowski

I think (but I am no expert) your first use of the metric (N = cross(G1, G2) ) is correct but its giving you the normal of the cell and depends on the orientation of the nodes.

If you just want a Function with the outward facing normal at the centre of each element for that particular mesh (but without altering the orientation) then you could create it as below.

from dolfin import *
from ufl import Jacobian
import numpy,dolfin as df
mesh = BoxMesh(Point(0, 0, 0), Point(2, 2, 2), 2, 2, 2)
mesh.coordinates()[:,0] *= .5   
mesh.coordinates()[:,1] *= .5
srf = BoundaryMesh(mesh, "exterior", order=True) 
J = Jacobian(srf)
TangentBasis1 = J[:,0]
TangentBasis2 = J[:,1]
N = cross(TangentBasis1, TangentBasis2) #  scaled cell normal field
Nh=project(N, VectorFunctionSpace(srf, 'DG', degree=0, dim=3))
Nv=numpy.array(Nh.vector()) 

V= df.TensorFunctionSpace(srf, 'DG', 0, shape=(3,))
fNormal =df.Function(V)
ret_dofmap = V.dofmap()
temparray=numpy.zeros((len([c for c in cells(srf)]),3)).flatten()
meshcentre=numpy.average(srf.coordinates(),axis=0)
print('mesh centre',meshcentre)
for c,mesh_cell in enumerate(cells(srf)):
    normal=numpy.array([v for v in mesh_cell.cell_normal()])
    normal2= normal*(1 if normal.dot(numpy.array([v for v in mesh_cell.midpoint()])-meshcentre)>0 else -1)
    print(c, 'centre',[v for v in mesh_cell.midpoint()],'actual cell normal',normal,'actual cell normal from metric',Nv[3*c:3*c+3],'corrected cell normal',normal2)
    temparray  [ret_dofmap.cell_dofs(mesh_cell.index())] =normal2.flatten()

fNormal.vector()[:] =temparray
deg=0
XDMFFile('norm2.deg.'+str(deg)+'.xdmf').write(project(fNormal, VectorFunctionSpace(srf, 'DG', degree=deg, dim=3)))

If you want to reorient the cells of the mesh and can’t work out how to use the fenics init_cell_orientations() function then almost any FE package should be able to make a mesh and orient the normals as you wish (you would then need to import the mesh).

One other thing I noticed - you project to DG2 to plot the normals , I thought you might be better projecting the normal field to DG 0 ? I find plotting centre of element vectors in paraview can be very confusing because the vectors at the centres of the elements seem to get extrapolated and averaged out to the nodes, I am not sure if its paraview or dolfin which does that.

3

Hi @Thomas2 ,
Thanks for your reply and suggestion. However, what I need is the (analytical) normal not at the center of each element, but at every dof (or integration point). I am solving some thin structure mechanics problems and have a parametric geometry defined on a UnitSquareMesh, thus the normal vector is not constant over an element. Would you happen to have any suggestions on how to achieve that with your type of approach?

Thanks again

4

HI Alexander,

… need is the (analytical) normal not at the center of each element, but at every dof (or integration point)…

Warning - I think I don’t fully understand what you need so my advice below might not be helpful;

I did try producing a normal on the function space using df.UserExpression, we actaully used DG0 - so only one node at the centre but I guess it would also work with (eg) CG1.

For us this approach works as long as the mesh was oriented so that the variation of the normals is everywhere smooth. But when we tried non manifolds with intersecting plates we couldn’t work out a robust way of getting this to work at the intersections

This gives you the average normal at each node rather than at integration points - I don’t really understand why you need the normals at each integration point so its probably not what you need.

class Metric(df.UserExpression):
    def __init__(self, mesh,sourceBasis=None,targetBasis=None,**kwargs):
        super(Metric, self).__init__( **kwargs)
        self.mesh = mesh
        self.targetBasis=targetBasis
        self.sourceBasis=sourceBasis

    def eval_cell(self, values, x, ufl_cell):
        mesh_cell = df.Cell(self.mesh, ufl_cell.index)
        xcell = numpy.asarray(mesh_cell.get_vertex_coordinates()).reshape(3, 3)
        e3=numpy.array([v for v in mesh_cell.cell_normal()])
        e3/=numpy.sqrt(numpy.dot(e3,e3))               
        e1=numpy.array( xcell[1]-xcell[0])
        e1/=numpy.sqrt(numpy.dot(e1,e1))  
        e2=numpy.cross(e1,e3)
        e2/=numpy.sqrt(numpy.dot(e2,e2))  
        basis=numpy.array([e1,e2,e3]).T
        mat=numpy.array([ [numpy.dot(basis[:,ii],self.sourceBasis[:,jj] ) for jj in range(3)] for ii in range(3)])
        values[:]=(-mat).flatten()

    def value_shape(self): return (3, 3)

# Tried both interpolate and project - results were similar with both of them.
#metric_at_nodes=df.Function( df.TensorFunctionSpace(mesh, 'DG', 0, shape=(3,3)))
#metric_at_nodes.interpolate(Metric(mesh,sourceBasis=numpy.identity(3)))    

metric_at_nodes=df.project(Metric(mesh,sourceBasis=numpy.identity(3)),V=df.TensorFunctionSpace(mesh, 'DG', 0, shape=(3,3)))

Because it doesn’t work with non manifold meshes we moved to the approach in my reply to your other question… that approach seems to work for intersecting plates - we integrate the strain energy over each element then use F=df.derivative(totalEnergy, uFunction, uTestFunction) - but we are still testing so not sure yet if it works for all scenarios.

solving some thin structure mechanics problems

Good to know you are using Fenics for thin structural mechanics - that it also what we are using it for. So far we have adapted Jack Hale’s “fenics shells” 2D python code for an Arnold-Brezzi partial reduced integration approach of a nonlinear Naghdi-Ressiner shell to 3D and added a Lagrange multiplier to link the drilling dof to the in plane shear strains. We aim to couple this to some 2nd order potential flow fluid code we have but its a stop start project so we probably won’t get to that part till next year or later. What are your thin structures?

1 Like
5

Hi Thomas,
Thanks for sharing these code snippets (especially the one from my other question - and yes, I too have found tIGAr very helpful). Indeed it seems I am attempting something similar. I’m working on structural pressurized membranes, currently investigating variational wrinkling models and trying to improve upon current techniques as part of my phd research. Previously I was only using geometrically exact models (like in Fenics Shells), but now I’d like to also consider arbitrary geometries, which necessitates making my own custom manifold meshes. We’re also interested in coupling these models with some flow solvers. I previously spent some time working to implement ALE methods in Tormod Landet’s Ocellaris to my problem, but I wasn’t making much progress. I plan to revisit that soon.