[petsc-users] How to construct DMPlex of cells with different topological dimension?

Matthew Knepley knepley at gmail.com
Wed Nov 3 06:19:36 CDT 2021


On Tue, Nov 2, 2021 at 10:49 PM 袁煕 <yuanxi at advancesoft.jp> wrote:

> Well! All objects in this world could definitely be modeled as 3D volumes.
> But considering costs of meshing and following calculations, in most cases
> it is neither practical nor necessary. And in some specific cases, e.g., to
> model cell membranes, 2D membrane elements may behave better than 3D volume
> ones. It is therefore, in many cases, modeling should be much simplified
> and in some cases they could be modeled as 2D or 1D ones.
>

I have _never_ suggested using 3D elements, as you can see from my reply to
Nicolas. I said you can use _exactly_ the element you want. I am not sure
how to be more clear.

   Matt


> Following paper is obtained after just google "parachute, FSI"
> *A parallel 3D computational method for fluid-structure interactions in
> parachute systems
> <https://www.researchgate.net/publication/228719684_A_parallel_3D_computational_method_for_fluid-structure_interactions_in_parachute_systems>*
>
> You can see they use cable (1D) and membrane (2D) elements to model
> parachutes. Those public software, such as ABAQUS, ANSYS, NASTRAN, etc all
> allow the mixed usage of  topologically different elements. They tell the
> same story. It is the reality and I don't think such reality would change
> in the near future.
>
> Best regards
>
> Yuan
>
> 2021年11月2日(火) 18:22 Matthew Knepley <knepley at gmail.com>:
>
>> On Mon, Nov 1, 2021 at 11:20 PM 袁煕 <yuanxi at advancesoft.jp> wrote:
>>
>>> Dear Matthew,
>>>
>>> I built a test problem using the strategy you suggested. It works! It is
>>> enough for me right now. Thank you very much.
>>>
>>> !       9----------8----------13
>>> !      /|            /|            /|
>>> !     / |           / |           / |
>>> !    /  |          /  |          /  |
>>> !   6---------7---------12  |
>>> !   |   |         |   |         |   |
>>> !   |   |         |   |         |   |
>>> !   |   |         |   |         |   |
>>> !   |   |         |   |         |   |
>>> !   |   5------|---4-------|--11--------17------16
>>> !   |  /          |  /          |  /            /           /
>>> !   | /           | /           | /            /           /
>>> !   |/            |/            |/            /           /
>>> !   2-- ------3---------10--------14-------15
>>>
>>>       coneSize = (/ 8,8,8,8, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 /)
>>>       cones = (/ 2,5,4,3,6,7,8,9,  3,4,11,10,7,12,13,8,  &
>>>            10,14,17,11, 10,14,17,11, 14,15,16,17, 14,15,16,17 /)
>>>
>>> There is still a problem left, however, that is how to build a 3D cell
>>> from a 2 vertex segment. I think I could construct a virtual tetrahedron.
>>>
>>> Finally, It can't be denied that this approach needs additional effort
>>> in mesh construction (and It is something strange to construct a 3D object
>>> from a one dimensional segment). In fact, structures with different
>>> topological dimensions are not that rare. Parachute, for example, may be
>>> composed with two dimensional parafoil and one dimensional cord. FRP(Fiber
>>> reinforced plastics) may be modelled by one dimensional reinforcement and
>>> three dimensional plastics. It is therefore, I am wondering, if PETSc would
>>> take into consider of such cases by, for example
>>>
>>> -   Enable a face, DM_POLYTOPE_SEGMENT2D, defined by one edge. And
>>> -   Eable cells, such
>>> as DM_POLYTOPE_QUADRILATERAL3D, DM_POLYTOPE_TRIANGLE3D and
>>> DM_POLYTOPE_SEGMENT3D, with one face.
>>>
>>
>> I think I am still not being completely clear. These types of cells
>> definitely exist. However, let's take the case of the parachute. The one
>> dimensional cords are part of
>> the two dimensional patches, which are part of a three dimensional
>> volume. This is how we have modeled it. You assemble different equations on
>> different pieces, but
>> that does not affect the mesh.
>>
>>   Thanks,
>>
>>      Matt
>>
>>
>>> Sorry for a newbie in PETSc to provide such a suggestion. But if it
>>> could be accomplished, it would help mech structural engineering, civil
>>> engineering, and solid mechanics applicationers.
>>>
>>> Best regards,
>>>
>>> Yuan
>>>
>>>
>>>
>>> 2021年11月1日(月) 18:32 Matthew Knepley <knepley at gmail.com>:
>>>
>>>> On Sun, Oct 31, 2021 at 12:21 PM 袁煕 <yuanxi at advancesoft.jp> wrote:
>>>>
>>>>> Dear Matt
>>>>>
>>>>> Thank you for your detailed explanation.
>>>>>
>>>>> First, I would like to answer your question about my application where
>>>>> low dimensional features are not embedded in the larger volume. It is quite
>>>>> general in structural engineering. For example, buildings are generally
>>>>> modelled as shells and beams, which are two and one dimension respectively.
>>>>> While  building foundation is modelled by solid elements, which is
>>>>> three dimension, at the same time.
>>>>>
>>>>
>>>> I think I see what you want now. Let me make a suggestion (along the
>>>> lines of what Mark said), and attempt to justify it by answering some
>>>> questions.
>>>>
>>>> Suggestion: I think you should consider using a volumetric mesh for
>>>> your problem
>>>>
>>>> Q1: Can I get the same algebraic system this way?
>>>>
>>>> Yes. You would only assign unknowns to faces and edges of the mesh
>>>> where you have shell and beam elements.
>>>>
>>>> Q2: What are the advantages?
>>>>
>>>> You would get topological connectivity of all parts of the structure,
>>>> automatic coupling of the different formulations,
>>>> and you could separately solve the different structures using block
>>>> preconditioners, while still forming a unified
>>>> residual so that you can guarantee a consistent solution.
>>>>
>>>> Q3: Wouldn't this be expensive?
>>>>
>>>> No. For the shells and beams, you would still need the vertices, edges,
>>>> and faces. First, by the Euler relation, these would
>>>> outnumbers the additional cells. Second, since no unknowns would be
>>>> associated with the cells, only additional memory in
>>>> the mesh would be used, not the system solves. Memory and time are
>>>> dominated by the solve, so this should be in the noise.
>>>>
>>>>   What do you think?
>>>>
>>>>      Thanks,
>>>>
>>>>         Matt
>>>>
>>>>
>>>>> Secondly, It is regrettably to find that DMComposite is not available
>>>>> for me, because all my solid, shell, and beam elements are connected each
>>>>> other.
>>>>>
>>>>> At last, I have build a simple program to see if DMPlexSetCellType()
>>>>> works for me, following the suggestion of functions in PETSc like
>>>>> DMPlexCreateCGNS. But it failed when it tried to do DMPlexInterpolate
>>>>> !       9----------8---------13
>>>>> !      /|            /|            /|
>>>>> !     / |           / |           / |
>>>>> !    /  |          /  |          /  |
>>>>> !   6---------7---------12  |
>>>>> !   |   |         |   |         |   |
>>>>> !   |   |         |   |         |   |
>>>>> !   |   |         |   |         |   |
>>>>> !   |   |         |   |         |   |
>>>>> !   |   5------|---4-------|-11--------17--------16
>>>>> !   |  /          |  /          |  /            /           /
>>>>> !   | /           | /           | /            /           /
>>>>> !   |/            |/            |/            /           /
>>>>> !   2---------3---------10--------14-------15
>>>>>
>>>>> The calculation result are follows. It seems that the cell type are
>>>>> set correctly but depth is still 2.
>>>>> --------------------------------------------------------------------
>>>>> DM Object: TestMesh 2 MPI processes
>>>>>   type: plex
>>>>> TestMesh in 3 dimensions:
>>>>>   0-cells: 16 0
>>>>>   3-cells: 20 (18) 0
>>>>> Labels:
>>>>>   celltype: 3 strata with value/size (7 (2), 4 (2), 0 (16))
>>>>>   depth: 2 strata with value/size (0 (16), 1 (20))
>>>>> [0]PETSC ERROR: --------------------- Error Message
>>>>> --------------------------------------------------------------
>>>>> [0]PETSC ERROR: Object is in wrong state
>>>>> [0]PETSC ERROR: Array was not checked out
>>>>> [0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble
>>>>> shooting.
>>>>> [0]PETSC ERROR: Petsc Development GIT revision:
>>>>> v3.16.0-351-g743e004674  GIT Date: 2021-10-29 09:32:23 -0500
>>>>> [0]PETSC ERROR: ./ex3f90 on a arch-linux-c-debug named DESKTOP-9ITFSBM
>>>>> by hillyuan Mon Nov  1 00:26:39 2021
>>>>> [0]PETSC ERROR: Configure options --with-cc=mpiicc --with-cxx=mpiicpc
>>>>> --with-fc=mpiifort --with-fortran-bindings=1
>>>>> --with-blaslapack-dir=/opt/intel/oneapi/mkl/2021.3.0
>>>>> --with-mkl_pardiso-dir=/opt/intel/oneapi/mkl/2021.3.0
>>>>> [0]PETSC ERROR: #1 DMRestoreWorkArray() at
>>>>> /home/hillyuan/programs/petsc/src/dm/interface/dm.c:1580
>>>>> [0]PETSC ERROR: #2 DMPlexRestoreRawFaces_Internal() at
>>>>> /home/hillyuan/programs/petsc/src/dm/impls/plex/plexinterpolate.c:323
>>>>> [0]PETSC ERROR: #3 DMPlexInterpolateFaces_Internal() at
>>>>> /home/hillyuan/programs/petsc/src/dm/impls/plex/plexinterpolate.c:375
>>>>> [0]PETSC ERROR: #4 DMPlexInterpolate() at
>>>>> /home/hillyuan/programs/petsc/src/dm/impls/plex/plexinterpolate.c:1340
>>>>>
>>>>> -----------------------------------------------------------------------------------------
>>>>> I attached my test program in this mail. It is much appreciated that
>>>>> you could provide any suggestion.
>>>>>
>>>>> Best regards,
>>>>>
>>>>> Yuan
>>>>>
>>>>>
>>>>> 2021年10月31日(日) 21:16 Matthew Knepley <knepley at gmail.com>:
>>>>>
>>>>>> On Thu, Oct 28, 2021 at 10:48 PM 袁煕 <yuanxi at advancesoft.jp> wrote:
>>>>>>
>>>>>>> Dear Matt,
>>>>>>>
>>>>>>>     My mesh is something like the following figure, which is
>>>>>>> composed of three elements : one hexahedron(solid element), one
>>>>>>> quadrilateral (shell element), and one line (beam element). I found the
>>>>>>> function "TestEmptyStrata" in file \dm\impls\plex\tests\ex11.c
>>>>>>> would be a good example to read in such a kind of mesh by using
>>>>>>> DMPlexSetCone. But a problem is that you should declare all faces and edges
>>>>>>> of hexahedron element, all edges in quadrilateral  element by
>>>>>>> DMPlexSetCone, otherwise PETsc could not do
>>>>>>> topological interpolation afterwards. Am I right here?
>>>>>>>    As general in FEM mesh, my mesh does not contain any information
>>>>>>> about faces or edges of solid elements. That's why I consider using
>>>>>>> DMCOMPOSITE. That is
>>>>>>>
>>>>>>> -   Put hexahedron, quadrilateral, and line elements into
>>>>>>> different  DM structures.
>>>>>>> -   do topological interpolation in those DMs separately.
>>>>>>> -   composite them.
>>>>>>>
>>>>>>>    Is there anything wrong in my above consideration?  Any
>>>>>>> suggestions?
>>>>>>>
>>>>>>>         ------------
>>>>>>>       /|             /|
>>>>>>>      / |           /  |  cell 0: Hex
>>>>>>>     /  |          /   |
>>>>>>>    ------------/   |
>>>>>>>    |   |         |   |
>>>>>>>    |   |         |   |   cell 1: Quad
>>>>>>>    |   --------|---|------------
>>>>>>>    |  /          |  /             /
>>>>>>>    | /           | /             /
>>>>>>>    |/            |/             /
>>>>>>>    -------------------------------------------
>>>>>>>                                   cell 2: line
>>>>>>>
>>>>>>> Much thanks for your help.
>>>>>>>
>>>>>>
>>>>>> If you are solving something where everything is embedded in a
>>>>>> volumetric mesh, then there is no problem. However, if you really have
>>>>>> the mesh above, where lower dimensional pieces are sticking out of
>>>>>> the mesh, then Plex can represent the mesh, but automatic interpolation
>>>>>> (creation of edges and faces) will not work. Why is this? We use
>>>>>> depth in the DAG as a proxy for cell dimension, but this will no longer work
>>>>>> if faces are not part of a volume.
>>>>>>
>>>>>> Will DMCOMPOSITE do what you want? It depends. It will be able to lay
>>>>>> out a vector, but it will not know about any topological connectivity
>>>>>> between the meshes and will not preallocate a Jacobian with any
>>>>>> interaction. If the meshes are truly separate, this is fine. If not, it is
>>>>>> not that
>>>>>> useful.
>>>>>>
>>>>>> Could you modify the existing code to support this? Yes, it would not
>>>>>> be terribly difficult. When you load the mesh, you must know what kind
>>>>>> of cell you are loading. You could explicitly set this using
>>>>>> DMPlexSetCellType(). Then, instead of taking a certain height stratum of
>>>>>> the DAG
>>>>>> to loop over, you would instead use all cells marked with a certain
>>>>>> cell type. The rest of the interpolation code should work fine.
>>>>>>
>>>>>> What kind of physics do you have where low dimensional features are
>>>>>> not embedded in the larger volume?
>>>>>>
>>>>>>   Thanks,
>>>>>>
>>>>>>      Matt
>>>>>>
>>>>>>
>>>>>>> Yuan
>>>>>>>
>>>>>>> 2021年10月28日(木) 22:05 Matthew Knepley <knepley at gmail.com>:
>>>>>>>
>>>>>>>> On Thu, Oct 28, 2021 at 4:59 AM 袁煕 <yuanxi at advancesoft.jp> wrote:
>>>>>>>>
>>>>>>>>> Dear Matt,
>>>>>>>>>
>>>>>>>>> Thank you for your quick response.
>>>>>>>>>
>>>>>>>>> I think what you mean is to build DAG from my mesh at first and
>>>>>>>>> then call DMPlexCreateFromDAG
>>>>>>>>> <https://petsc.org/main/docs/manualpages/DMPLEX/DMPlexCreateFromDAG.html#DMPlexCreateFromDAG>()
>>>>>>>>> to construct DMPlex.
>>>>>>>>>
>>>>>>>>
>>>>>>>> No, I do not mean that.
>>>>>>>>
>>>>>>>>
>>>>>>>>> A new problem is,  as I know, the function DMPlexInterpolate would
>>>>>>>>> generate points with different depth. What's the difference  between those
>>>>>>>>> faces and segment elements generated by  DMPlexInterpolate  with that
>>>>>>>>> defined by the original mesh, or should we not use DMPlexInterpolate in
>>>>>>>>> such a case?
>>>>>>>>>
>>>>>>>>> On the other hand, can DMComposite be used in this case by
>>>>>>>>> defining DMPlex with different topological dimensions at first and then
>>>>>>>>> composite them?
>>>>>>>>>
>>>>>>>>
>>>>>>>> You do not need that. I am obviously not understanding your
>>>>>>>> question. My short answer is that Plex _already_ handles cells of different
>>>>>>>> dimension automatically without anything extra.
>>>>>>>>
>>>>>>>> Maybe it would help if you defined a specific problem you have.
>>>>>>>>
>>>>>>>>   Thanks,
>>>>>>>>
>>>>>>>>       Matt
>>>>>>>>
>>>>>>>>
>>>>>>>>> Thanks in advance.
>>>>>>>>>
>>>>>>>>> Yuan
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> 2021年10月27日(水) 19:27 Matthew Knepley <knepley at gmail.com>:
>>>>>>>>>
>>>>>>>>>> On Wed, Oct 27, 2021 at 4:50 AM 袁煕 <yuanxi at advancesoft.jp> wrote:
>>>>>>>>>>
>>>>>>>>>>> Hi,
>>>>>>>>>>>
>>>>>>>>>>> I am trying to parallelize my serial FEM program using PETSc.
>>>>>>>>>>> This program calculates structure deformation by using various types of
>>>>>>>>>>> elements such as solid, shell, beam, and truss. At the very beginning, I
>>>>>>>>>>> found it was hard for me to put such kinds of elements into DMPlex. Because
>>>>>>>>>>> solid elements are topologically three dimensional, shell element two, and
>>>>>>>>>>> beam or truss are topologically one-dimensional elements. After reading
>>>>>>>>>>> chapter 2.10: "DMPlex: Unstructured Grids in PETSc" of users manual
>>>>>>>>>>> carefully,  I found the provided functions, such as DMPlexSetCone, cannot
>>>>>>>>>>> declare those topological differences.
>>>>>>>>>>>
>>>>>>>>>>> My question is : Is it possible and how to define all those
>>>>>>>>>>> topologically different elements into a DMPlex struct?
>>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> Yes. The idea is to program in a dimension-independent way, so
>>>>>>>>>> that the code can handle cells of any dimension.
>>>>>>>>>> What you probably want is the "depth" in the DAG representation,
>>>>>>>>>> which you can think of as the dimension of a cell.
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> https://petsc.org/main/docs/manualpages/DMPLEX/DMPlexGetPointDepth.html#DMPlexGetPointDepth
>>>>>>>>>>
>>>>>>>>>>   Thanks,
>>>>>>>>>>
>>>>>>>>>>      Matt
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>> Thanks in advance!
>>>>>>>>>>>
>>>>>>>>>>> Best regards,
>>>>>>>>>>>
>>>>>>>>>>> Yuan.
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> --
>>>>>>>>>> What most experimenters take for granted before they begin their
>>>>>>>>>> experiments is infinitely more interesting than any results to which their
>>>>>>>>>> experiments lead.
>>>>>>>>>> -- Norbert Wiener
>>>>>>>>>>
>>>>>>>>>> https://www.cse.buffalo.edu/~knepley/
>>>>>>>>>> <http://www.cse.buffalo.edu/~knepley/>
>>>>>>>>>>
>>>>>>>>>
>>>>>>>>
>>>>>>>> --
>>>>>>>> What most experimenters take for granted before they begin their
>>>>>>>> experiments is infinitely more interesting than any results to which their
>>>>>>>> experiments lead.
>>>>>>>> -- Norbert Wiener
>>>>>>>>
>>>>>>>> https://www.cse.buffalo.edu/~knepley/
>>>>>>>> <http://www.cse.buffalo.edu/~knepley/>
>>>>>>>>
>>>>>>>
>>>>>>
>>>>>> --
>>>>>> What most experimenters take for granted before they begin their
>>>>>> experiments is infinitely more interesting than any results to which their
>>>>>> experiments lead.
>>>>>> -- Norbert Wiener
>>>>>>
>>>>>> https://www.cse.buffalo.edu/~knepley/
>>>>>> <http://www.cse.buffalo.edu/~knepley/>
>>>>>>
>>>>>
>>>>
>>>> --
>>>> What most experimenters take for granted before they begin their
>>>> experiments is infinitely more interesting than any results to which their
>>>> experiments lead.
>>>> -- Norbert Wiener
>>>>
>>>> https://www.cse.buffalo.edu/~knepley/
>>>> <http://www.cse.buffalo.edu/~knepley/>
>>>>
>>>
>>
>> --
>> What most experimenters take for granted before they begin their
>> experiments is infinitely more interesting than any results to which their
>> experiments lead.
>> -- Norbert Wiener
>>
>> https://www.cse.buffalo.edu/~knepley/
>> <http://www.cse.buffalo.edu/~knepley/>
>>
>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20211103/cf544492/attachment-0001.html>


More information about the petsc-users mailing list