Hello All,
I’m working on a model of a composite under pressure, but have run into some difficulties. My model is based on the hyper elasticity tutorial, but has additional terms to model the compaction of the fibres within the material. The fibre lay is controlled by the vectors e0 and n0 (where e0 is the direction of the fibre and n0 is normal to it). The model will run with a simple block, but in order to use more complex geometry, I need to compute e0 and n0 for each cell, so the fibres follow the curve of the material (see image for example).
Any help or pointers would be greatly appreciated!
There is a quite similar question here. The simplest solution would be, if you know the directions e^0 and n^0 analytically, to just specify them in terms of spatial coordinates. For instance, if they are (normalized) 2D polar basis vectors (as suggested by the figure), and the origin of your coordinate system is the center of the circle, you could do something like:
# Cartesian coordinates
x = SpatialCoordinate(mesh)
# Unit vector in radial direction
n0 = x/sqrt(dot(x,x))
# 90-degree rotation
e0 = as_vector([n0[1],-n0[0]])
In general, I would argue that trying to determine material directions from the mesh is not best practice, since they are part of the mathematical problem statement, not the discretization, and should be well-defined prior to generating a mesh.
3 Likes
My apologies for not replying sooner, but thank you so much for your answer. I’ve got this to run successfully in my code, but it does make for a computationally-expensive model. Do you have any suggestions for how I might reduce the run-time of this model?
Thanks again.
It really depends on where (assembly, solution, etc.) the time is being spent. I can think of one potential issue related to this discussion, which would affect system assembly time. Forms with complicated integrands (e.g., involving transformations to different coordinate systems, as above) often cause FEniCS’s default estimation of quadrature degree to blow up, leading to excessive numbers of quadrature points. You can manually set quadrature degree with
dx = dx(metadata={"quadrature_degree":n})
(and likewise for ds and dS), where n is an integer giving a custom polynomial degree up to which the quadrature is exact.
1 Like
Thanks, I’ll try adding this into my code!
Thank you again for all your help, @kamensky!
