adding mesh elements over time

Asked by Nick Davies

Hi everyone,
I am trying to figure out a way to have my mesh grow in time. Basically what it is doing is the original mesh is being deformed (compressed) then a new layer is added with some different values such as Youngs modulus etc to the outside of the compressed mesh, however the new elements on the outside are not under compression. After this the new mesh which includes the original compressed mesh and the added uncompressed elements are deformed again and then the next layer is added and so on. It is to simulate the growth of a tree stem or branch.
What this is emulating in the physical world is as a tree grows the cells near the bottom of the stem are compressed by the trees self weight. Because the tree is now shorter as it produces new cells on the outside of the compressed cells they are not compressed as they have little weight above them. What this results in is the cells in the center being further compressed than the cells on the outer surface of the stem because the cells are younger at the out edge.

I hope this makes some sense to everyone. Does anybody have any ideas as to how to implement this?

Thanks for any help

Cheers Nick

P.S. I am using the dev of fenics and generating the mesh with the cone() command

Question information

Language:
English Edit question
Status:
Solved
For:
FEniCS Project Edit question
Assignee:
No assignee Edit question
Solved by:
Nick Davies
Solved:
Last query:
Last reply:
Revision history for this message
Johan Hake (johan-hake) said :
#1

Adding or removing (in place) cells from already existing meshes is not
possible at the moment.

You can however create a new mesh with the added cells. You then need to
add logic how your transfer functional information between the two pdes.
This is what the whole adapt interface of dolfin does.

If you want to try this approach, you need to figure out what vertices
and cells need to be added. Open a new empty mesh in MeshEditor, and
initialize it with the updated num_vertices and num_cells. Then you add
information about the old mesh, together with the new vertices and cells.

You should be able to handle the resulting mesh hierarchy with the
Hierarchical interface of Mesh, ie. set_parent, and set_child.

Johan

On 12/03/2012 12:55 AM, Nick Davies wrote:
> New question #215827 on FEniCS Project:
> https://answers.launchpad.net/fenics/+question/215827
>
> Hi everyone, I am trying to figure out a way to have my mesh grow in
> time. Basically what it is doing is the original mesh is being
> deformed (compressed) then a new layer is added with some different
> values such as Youngs modulus etc to the outside of the compressed
> mesh, however the new elements on the outside are not under
> compression. After this the new mesh which includes the original
> compressed mesh and the added uncompressed elements are deformed
> again and then the next layer is added and so on. It is to simulate
> the growth of a tree stem or branch. What this is emulating in the
> physical world is as a tree grows the cells near the bottom of the
> stem are compressed by the trees self weight. Because the tree is now
> shorter as it produces new cells on the outside of the compressed
> cells they are not compressed as they have little weight above them.
> What this results in is the cells in the center being further
> compressed than the cells on the outer surface of the stem because
> the cells are younger at the out edge.
>
> I hope this makes some sense to everyone. Does anybody have any ideas
> as to how to implement this?
>
> Thanks for any help
>
> Cheers Nick
>
> P.S. I am using the dev of fenics and generating the mesh with the
> cone() command
>

Revision history for this message
Nick Davies (ntd14) said :
#2

Thanks for that Johan,
I have been tinkering with the approach you described above over the last couple days and an hoping someone can help with a problem which is arising. I am using the following code segment to retrieve the solution to the displacement which is calculated using solve() and converting it into a numpy array.

solve(F == 0, u, bc, J=J)
pp = u.vector()
u_array = pp.array()

I have used various methods to try to use this array to manipulate the mesh coordinates generated by the mesh.coordinates() command. From a visual inspection it appears for the most part the u_array is structured as z,y,x,z,y,x,...... however this does not always seam to be the case, with random points throughout the vector appearing to be out of order. This is because I would expect the displacements in the horizontal directions to be small and the vertical displacements to be large as the compression is currently set to use a Poisson ratio of 0. For about 95% of the entries in u_array the z,y,x patten follows this expectation with the vertical (z) displacements being ~1000 times grater than the horizontal (x,y) displacements but the rest of the entries seam to disrupt this with some entries in the z direction being very small (same order as the x and y entries) and some x and y entries being large (same order as the z entries).

Is this a good method for extracting the values from u or is there are more sensible way of doing it?

Thanks for any help

Revision history for this message
Johan Hake (johan-hake) said :
#3

On 12/06/2012 02:11 AM, Nick Davies wrote:
> Question #215827 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/215827
>
> Nick Davies posted a new comment:
> Thanks for that Johan,
> I have been tinkering with the approach you described above over the last couple days and an hoping someone can help with a problem which is arising. I am using the following code segment to retrieve the solution to the displacement which is calculated using solve() and converting it into a numpy array.
>
> solve(F == 0, u, bc, J=J)
> pp = u.vector()
> u_array = pp.array()
>
> I have used various methods to try to use this array to manipulate the
> mesh coordinates generated by the mesh.coordinates() command. From a
> visual inspection it appears for the most part the u_array is structured
> as z,y,x,z,y,x,...... however this does not always seam to be the case,
> with random points throughout the vector appearing to be out of order.

That is correct. You cannot expect anything about the ordering. You can
however build a map between coordinate. I have made such a function for
CG1 space and CG1 vector space:

def dofs_to_verts(mesh):
    """
    Compute the dof to vert mapping
    """
    import dolfin as d

    coords = mesh.coordinates()

    # FunctionSpaces
    V = d.FunctionSpace(mesh, "CG", 1)
    Vv = d.VectorFunctionSpace(mesh, "CG", 1)

    # Dof arrays
    dofs_to_vert = np.zeros(mesh.num_vertices(), dtype=np.uintp)
    vectordofs_to_vert = np.zeros((mesh.num_vertices()*3), dtype=np.uintp)
    vectordofs_to_subvert = np.zeros((mesh.num_vertices()*3),
dtype=np.uintp)

    # DofMaps
    dm = V.dofmap()
    dms = [Vv.sub(i).dofmap() for i in range(3)]

    # Iterate over cells and collect dofs
    for cell in cells(mesh):
        cell_ind = cell.index()
        vert_inds = cell.entities(0)

        # Check that dof coordinates are the same as the vertex coordinates
        assert(np.all(coords[vert_inds]==dm.tabulate_coordinates(cell)))
        dofs_to_vert[dm.cell_dofs(cell_ind)] = vert_inds

        # Iterate over the sub dof maps
        for i, (dms_i, dmcs_i) in enumerate(zip(dms, dmcs)):
            vectordofs_to_vert[dms_i.cell_dofs(cell_ind)] = vert_inds
            vectordofs_to_subvert[dms_i.cell_dofs(cell_ind)] = i

    # Return dof to vertex mapping
    return dofs_to_vert, vectordofs_to_vert, vectordofs_to_subvert

Feel free to use and abuse. Note it is hard coded for 3D meshes.

> This is because I would expect the displacements in the horizontal
> directions to be small and the vertical displacements to be large as the
> compression is currently set to use a Poisson ratio of 0. For about 95%
> of the entries in u_array the z,y,x patten follows this expectation
> with the vertical (z) displacements being ~1000 times grater than the
> horizontal (x,y) displacements but the rest of the entries seam to
> disrupt this with some entries in the z direction being very small (same
> order as the x and y entries) and some x and y entries being large (same
> order as the z entries).

I cannot say anything about the expected outcome of your PDE as no code
is provided. I am not that sturdy in mechanics, but others are.

> Is this a good method for extracting the values from u or is there are
> more sensible way of doing it?

Try the index mapping above. You most probably need the reversed mapping
though. From vertex to dof. But having the above mapping you should
easily get the reversed.

Johan

> Thanks for any help
>

Revision history for this message
Nick Davies (ntd14) said :
#4

Thanks for the code it helped a lot.