[petsc-users] Additiional DoF per cell

Fri May 8 07:01:43 CDT 2026

Moving on to the values of the closure.

Initially I create a vector for the DM with DMCreateLocal/GlobalVector, which is used for e.g. the solution. Following the example before, it will have size 8 (coordinates) + nc (for the p field).

I am getting a seemingly random order when I test with different values of nc. Say, I initialize the vector (VecSetValues) with the coordinates x_i, y_i . The rest nc components are set to zero. The closure returns:

-- nc = 1
closure = y0 x1 y1 x2 y2 x3 y3 0 x0

-- nc = 3
closure = y1 x2 y2 x3 y3 0 0 0 x0 y0 x1

-- nc = 7
closure = y3 0 0 0 0 0 0 0 x0 y0 x1 y1 x2 y2

Looking at that output:
- The first value is always a y coordinate (second component of u field).
- The position of the nc = 0 components change (starting with nc = 8, all zero values come first, followed by the coordinates).

Looking again at ex77.c (https://urldefense.us/v3/__https://petsc.org/release/src/snes/tutorials/ex77.c.html__;!!G_uCfscf7eWS!ZJz_pZR7M_mMi6gIDNMkNQy_6DO0K3sPQyCW-lCG3GWkYxEDauwvsaPpgiykvB00ywn3uvdURZsfZd4cNU2ZdMtLIvmj8bxC$ ), it says the closure stacks up the fields, which aligns with the output above (all u values, then all p values), but which is the starting point seems arbitrary.

The vector is initialized as [x0 y0 .... x3 y3 0 0 ... 0], which may not be correct to start with.

I noticed I could instead use DMPlexVecGetClosureAtDepth( ..., depth = 0, ...) so that I only get the u field DoF. But I will need to deal with the p field values as well at some point.

Thanks,
Noam

On Friday, May 8th, 2026 at 11:05 AM, Matthew Knepley <knepley at gmail.com> wrote:

> On Fri, May 8, 2026 at 5:08 AM Noam T. <dontbugthedevs at proton.me> wrote:
>
>> On the plus side, using the flags you provided in the previous email
>>
>> -dm_plex_simplex 0 -mech_petscspace_degree 1 -pres_petscspace_degree 1 -pres_petscspace_poly_tensor 0 -pres_petscdualspace_lagrange_continuity 0
>>
>> the output of -ds_view is identical to yours, so setting up the fields should be correct.
>> Thanks. I forgot about the flag "-[]_petscdualspace_lagrange_continuity 0" so as to have a discontinuous field.
>>
>> However, the closure size was still wrong.
>>
>> But, I figured out the reason: the manual set up of the section.
>> For the argument "numDof" in DMPlexCreateSection(), I was computing the number of DoF with PetscSectionGetDof() (which was fine until now, since there was a single field), instead of PetscSectionGetFieldDof().
>> So as you previously suggested
>>
>>> [...] It sounds like somehow you have a copy of your fields.
>>
>> the number of DoF per field were indeed "doubled" with the sum of both fields.
>
> Excellent. Let me know if you have any more issues.
>
> Thanks,
>
> Matt
>
>> Thanks,
>> Noam
>>
>> On Wednesday, May 6th, 2026 at 2:06 PM, Matthew Knepley <knepley at gmail.com> wrote:
>>
>>> Okay, this is a problem with specifying the space I believe. You want a discontinuous space for pressure.
>>>
>>> Here is what I get
>>>
>>> Discrete System with 2 fields
>>> cell total dim 11 total comp 3
>>> Field Q1 FEM 2 components (implicit) (Nq 4 Nqc 1) 1-jet
>>> PetscFE Object: Q1 (mech_) 1 MPI process
>>> type: vector
>>> Vector Finite Element in 2 dimensions with 2 components
>>> PetscSpace Object: Q1 (mech_) 1 MPI process
>>> type: sum
>>> Space in 2 variables with 2 components, size 8
>>> Sum space of 2 concatenated subspaces (all identical)
>>> PetscSpace Object: Q1 (mech_sumcomp_) 1 MPI process
>>> type: tensor
>>> Space in 2 variables with 1 components, size 4
>>> Tensor space of 2 subspaces (all identical)
>>> PetscSpace Object: sum component tensor component (mech_sumcomp_tensorcomp_) 1 MPI process
>>> type: poly
>>> Space in 1 variables with 1 components, size 2
>>> Polynomial space of degree 1
>>> PetscDualSpace Object: Q1 (mech_) 1 MPI process
>>> type: sum
>>> Dual space with 2 components, size 8
>>> Sum dual space of 2 concatenated subspaces (all identical)
>>> PetscDualSpace Object: Q1 1 MPI process
>>> type: lagrange
>>> Dual space with 1 components, size 4
>>> Continuous tensor Lagrange dual space
>>> Quadrature on a quadrilateral of order 3 on 4 points (dim 2)
>>> Field Q1 FEM 1 component (implicit) (Nq 4 Nqc 1) 1-jet
>>> PetscFE Object: Q1 (pres_) 1 MPI process
>>> type: basic
>>> Basic Finite Element in 2 dimensions with 1 components
>>> PetscSpace Object: Q1 (pres_) 1 MPI process
>>> type: poly
>>> Space in 2 variables with 1 components, size 3
>>> Polynomial space of degree 1
>>> PetscDualSpace Object: Q1 (pres_) 1 MPI process
>>> type: lagrange
>>> Dual space with 1 components, size 3
>>> Discontinuous Lagrange dual space
>>> Quadrature on a quadrilateral of order 3 on 4 points (dim 2)
>>> Weak Form System with 2 fields
>>>
>>> I have attached the code I used.
>>>
>>> Thanks,
>>>
>>> Matt
>>>
>>> On Wed, May 6, 2026 at 9:32 AM Noam T. <dontbugthedevs at proton.me> wrote:
>>>
>>>>>> Field p FEM 3 components (implicit) (Nq 4 Nqc 1) 1-jet
>>>>>
>>>>> 1. Why does pressure have 3 components in 2D?
>>>>
>>>> I used nc = 3 when calling PetscFECreateDefault() for the pressure field. Confusion between wanting 3 addditional DoF and components.
>>>>
>>>> With nc = 1 the closure size is 18 instead.
>>>>
>>>>> 2. Did you look at the incompressible example SNES ex69? It sets up this exact element.
>>>>
>>>> I changed some of my code following the example and got rid of some errors. But I'll have a more thorough look again.
>>>>
>>>>> 3. Please send me that code that sets up your DM. It sounds like somehow you have a copy of your fields.
>>>>>
>>>>> Thanks,
>>>>
>>>> For reference, -ds_view with only one (displacements) fields shows:
>>>>
>>>> ---
>>>> Discrete System with 1 fields
>>>> cell total dim 8 total comp 2
>>>> Field u FEM 2 components (implicit) (Nq 4 Nqc 1) 1-jet
>>>> PetscFE Object: u (disp_) 1 MPI process
>>>> type: vector
>>>> Vector Finite Element in 2 dimensions with 2 components
>>>> PetscSpace Object: Q1 (disp_) 1 MPI process
>>>> type: sum
>>>> Space in 2 variables with 2 components, size 8
>>>> Sum space of 2 concatenated subspaces (all identical)
>>>> PetscSpace Object: Q1 (disp_sumcomp_) 1 MPI process
>>>> type: tensor
>>>> Space in 2 variables with 1 components, size 4
>>>> Tensor space of 2 subspaces (all identical)
>>>> PetscSpace Object: sum component tensor component (disp_sumcomp_tensorcomp_) 1 MPI process
>>>> type: poly
>>>> Space in 1 variables with 1 components, size 2
>>>> Polynomial space of degree 1
>>>> PetscDualSpace Object: Q1 (disp_) 1 MPI process
>>>> type: sum
>>>> Dual space with 2 components, size 8
>>>> Sum dual space of 2 concatenated subspaces (all identical)
>>>> PetscDualSpace Object: Q1 1 MPI process
>>>> type: lagrange
>>>> Dual space with 1 components, size 4
>>>> Continuous tensor Lagrange dual space
>>>> Quadrature on a quadrilateral of order 3 on 4 points (dim 2)
>>>> Weak Form System with 1 fields---
>>>>
>>>> and -dm_view
>>>>
>>>> ---
>>>> DM Object: 1 MPI process
>>>> type: plex
>>>> DM_0x55c117db1b00_2 in 2 dimensions:
>>>> Number of 0-cells per rank: 4
>>>> Number of 1-cells per rank: 4
>>>> Number of 2-cells per rank: 1
>>>> Labels:
>>>> celltype: 3 strata with value/size (0 (4), 1 (4), 4 (1))
>>>> depth: 3 strata with value/size (0 (4), 1 (4), 2 (1))
>>>> Cell Sets: 1 strata with value/size (1 (1))
>>>> surf: 1 strata with value/size (1 (1))
>>>> Face Sets: 2 strata with value/size (2 (1), 3 (1))
>>>> edge_right: 1 strata with value/size (3 (1))
>>>> edge_left: 1 strata with value/size (2 (1))
>>>> Vertex Sets: 4 strata with value/size (4 (1), 5 (1), 6 (1), 7 (1))
>>>> P1: 1 strata with value/size (4 (1))
>>>> P2: 1 strata with value/size (5 (1))
>>>> P3: 1 strata with value/size (6 (1))
>>>> P4: 1 strata with value/size (7 (1))
>>>> Field u:
>>>> adjacency FEM
>>>> /// With the second field, the two lines below are added
>>>> Field p:
>>>> adjacency FEM
>>>> ---
>>>>
>>>> DM set up
>>>>
>>>> DM :: dm_mesh
>>>> PetscFE :: u_FE, p_FE
>>>> PetscInt :: dim, dim_u_FE, tdim
>>>> PetscDS :: dm_ds
>>>>
>>>> DMCreate(PETSC_COMM_WORLD, dm_mesh)
>>>> DMSetType(dm_mesh, DMPLEX)
>>>> DMGetDimension(dm_mesh, dim) // dim = 2
>>>> DMSetFromOptions(dm_mesh) // this was called twice by mistake, now removed; no changes that I can see
>>>>
>>>> -- using flags
>>>>
>>>> -dm_plex_filename QUAD1.msh
>>>> -dm_plex_interpolate 1
>>>> -dm_plex_gmsh_use_generic
>>>> -dm_plex_gmsh_use_regions
>>>> -dm_plex_gmsh_multiple_tags
>>>> -dm_plex_gmsh_mark_vertices
>>>>
>>>> -- some operations with labels
>>>>
>>>> // u field
>>>> PetscFECreateDefault(PETSC_COMM_SELF, dim, dim, PETSC_FALSE, "mech_", PETSC_DETERMINE, u_FE)
>>>> PetscFEGetDimension(u_FE, dim_u_FE) // dim_u_FE = 8
>>>> DMAddField(dm_mesh, NULL, (PetscObject)u_FE)
>>>>
>>>> // p field
>>>> PetscFECreateDefault(PETSC_COMM_SELF, dim , 1, PETSC_FALSE, "pres_", PETSC_DETERMINE, p_FE)
>>>> PetscFEGetDimension(p_FE, dim_p_FE) // dim_p_FE = 1
>>>> DMAddField(dm_mesh, NULL, (PetscObject)p_FE)
>>>>
>>>> DMCreateDS(dm_mesh)
>>>> DMGetDS(dm_mesh, dm_DS)
>>>> PetscDSGetTotalDimension(dm_DS, tdim) // tdim
>>>>
>>>> with flags
>>>>
>>>> -disp_petscdualspace_lagrange_node_type equispaced
>>>> -disp_petscdualspace_lagrange_node_endpoints 1
>>>>
>>>> I believe these are all the DM related calls.
>>>>
>>>> Thanks,
>>>> Noam
>>>>
>>>> On Wednesday, May 6th, 2026 at 11:39 AM, Matthew Knepley <knepley at gmail.com> wrote:
>>>>
>>>>> On Wed, May 6, 2026 at 6:56 AM Noam T. <dontbugthedevs at proton.me> wrote:
>>>>>
>>>>>> Hello,
>>>>>>
>>>>>>> That is wrong. Are you sure about 22?
>>>>>>
>>>>>> That's what I get, from the argument "csize" in DMPlexGetVecClosure().
>>>>>>
>>>>>> Here's the output of PetscDSView:
>>>>>>
>>>>>> Discrete System with 2 fields
>>>>>> cell total dim 11 total comp 5
>>>>>> Field u FEM 2 components (implicit) (Nq 4 Nqc 1) 1-jet
>>>>>> PetscFE Object: u (mech_) 1 MPI process
>>>>>> type: vector
>>>>>> Vector Finite Element in 2 dimensions with 2 components
>>>>>> PetscSpace Object: Q1 (mech_) 1 MPI process
>>>>>> type: sum
>>>>>> Space in 2 variables with 2 components, size 8
>>>>>> Sum space of 2 concatenated subspaces (all identical)
>>>>>> PetscSpace Object: Q1 (mech_sumcomp_) 1 MPI process
>>>>>> type: tensor
>>>>>> Space in 2 variables with 1 components, size 4
>>>>>> Tensor space of 2 subspaces (all identical)
>>>>>> PetscSpace Object: sum component tensor component (mech_sumcomp_tensorcomp_) 1 MPI process
>>>>>> type: poly
>>>>>> Space in 1 variables with 1 components, size 2
>>>>>> Polynomial space of degree 1
>>>>>> PetscDualSpace Object: Q1 (mech_) 1 MPI process
>>>>>> type: sum
>>>>>> Dual space with 2 components, size 8
>>>>>> Sum dual space of 2 concatenated subspaces (all identical)
>>>>>> PetscDualSpace Object: Q1 1 MPI process
>>>>>> type: lagrange
>>>>>> Dual space with 1 components, size 4
>>>>>> Continuous tensor Lagrange dual space
>>>>>> Quadrature on a quadrilateral of order 3 on 4 points (dim 2)
>>>>>> Field p FEM 3 components (implicit) (Nq 4 Nqc 1) 1-jet
>>>>>
>>>>> 1. Why does pressure have 3 components in 2D?
>>>>>
>>>>>> PetscFE Object: p (pres_) 1 MPI process
>>>>>> type: vector
>>>>>> Vector Finite Element in 2 dimensions with 3 components
>>>>>> PetscSpace Object: Q0 (pres_) 1 MPI process
>>>>>> type: sum
>>>>>> Space in 2 variables with 3 components, size 3
>>>>>> Sum space of 3 concatenated subspaces (all identical)
>>>>>> PetscSpace Object: Q0 (pres_sumcomp_) 1 MPI process
>>>>>> type: tensor
>>>>>> Space in 2 variables with 1 components, size 1
>>>>>> Tensor space of 2 subspaces (all identical)
>>>>>> PetscSpace Object: sum component tensor component (pres_sumcomp_tensorcomp_) 1 MPI process
>>>>>> type: poly
>>>>>> Space in 1 variables with 1 components, size 1
>>>>>> Polynomial space of degree 0
>>>>>> PetscDualSpace Object: Q0 (pres_) 1 MPI process
>>>>>> type: sum
>>>>>> Dual space with 3 components, size 3
>>>>>> Sum dual space of 3 concatenated subspaces (all identical)
>>>>>> PetscDualSpace Object: Q0 1 MPI process
>>>>>> type: lagrange
>>>>>> Dual space with 1 components, size 1
>>>>>> Discontinuous tensor Lagrange dual space
>>>>>> Quadrature on a quadrilateral of order 3 on 4 points (dim 2)
>>>>>> Weak Form System with 2 fields
>>>>>>
>>>>>> using as the second field (dim = 2, nc = 3, prefix = "pres") and setting the FE object name to "p".
>>>>>>
>>>>>> The entry "cell total dim 10" agrees with the output of PetscFEGetTotalDimension(); "comp = 5" I assume is 2 (u) + 3 (p)
>>>>>>
>>>>>> However, the closure of a Vec is still 22.
>>>>>
>>>>> 2. Did you look at the incompressible example SNES ex69? It sets up this exact element.
>>>>>
>>>>> 3. Please send me that code that sets up your DM. It sounds like somehow you have a copy of your fields.
>>>>>
>>>>> Thanks,
>>>>>
>>>>> Matt
>>>>>
>>>>>> On Tuesday, May 5th, 2026 at 1:58 PM, Matthew Knepley <knepley at gmail.com> wrote:
>>>>>>
>>>>>>> On Tue, May 5, 2026 at 9:06 AM Noam T. via petsc-users <petsc-users at mcs.anl.gov> wrote:
>>>>>>>
>>>>>>>> Hello,
>>>>>>>>
>>>>>>>> I am trying to work with a "mixed" FE discretization, where besides the usual displacements DoF in nodes, there is an additional field (e.g. pressure) that is also part of the system. This additional field has a certain number of additional dof : p0, p1, p2...
>>>>>>>>
>>>>>>>> Looking at example 77 ([https://urldefense.us/v3/__https://petsc.org/release/src/snes/tutorials/ex77.c.html*5D(https:/*urldefense.us/v3/__https:/*petsc.org/release/src/snes/tutorials/ex77.c.html__;!!G_uCfscf7eWS!aZNOxZZeoHYqe6lRZV0wHMVVTu3YIEfc1Sr-dca7xiGVfG3enULFD4_g6_fxWIa7A_Cu9LcXU-MuxRj4sOwDRyWqvsdGdrEd$)__;JS8v!!G_uCfscf7eWS!ZJz_pZR7M_mMi6gIDNMkNQy_6DO0K3sPQyCW-lCG3GWkYxEDauwvsaPpgiykvB00ywn3uvdURZsfZd4cNU2ZdMtLIqGUiP8y$ ) this seems to be handled with an additional field, added to the DM. I've tried so, but then I am getting from the DM arrays whose size/contents are not what I expected.
>>>>>>>>
>>>>>>>> For example, a four-noded Q1 element (PetscFE created with a 2x2 Gauss rule for quadrature), with just one field, a call to DMPlexGetVecClosure() gives me an array with 8 entries (say, the coordinates of the initial mesh): x = x0, y0, .... x3, y3
>>>>>>>>
>>>>>>>> Then add the new field:
>>>>>>>>
>>>>>>>> PetscFECreateDefault(..., nc = 3, ..., p_FE) /* not sure about the value of nc here */
>>>>>>>> DMAddField(dm, ..., p_FE)
>>>>>>>>
>>>>>>>> The total dimension, from PetscDSGetTotalDimension(), is now 11 (8 + 3). This results in a closure of size 22.
>>>>>>>
>>>>>>> That is wrong. Are you sure about 22?
>>>>>>>
>>>>>>>> However, what I am looking for is a closure of size 8 + 3 i.e. the original 8 DoF at the nodes, plus exactly 3 DoF (for the cell, so to speak) p0, p1, p2, so that in a system with "block" matrices of the form
>>>>>>>>
>>>>>>>> [K_uu, K_up | K_pu, K_pp] { u | p } = RHS
>>>>>>>>
>>>>>>>> the unknowns {p} has size 3 per cell. Is this possible? I tried some combinations of dim / nc for the new PetscFE, but when creating the DM section I get errors e.g.
>>>>>>>
>>>>>>> This is what you should get. I definitely have examples that do this. For example, here is Q1-P0 (I think that is what you are suggesting) for incompressible Stokes
>>>>>>>
>>>>>>> https://urldefense.us/v3/__https://gitlab.com/petsc/petsc/-/blob/main/src/snes/tutorials/ex69.c?ref_type=heads*L3454__;Iw!!G_uCfscf7eWS!ZJz_pZR7M_mMi6gIDNMkNQy_6DO0K3sPQyCW-lCG3GWkYxEDauwvsaPpgiykvB00ywn3uvdURZsfZd4cNU2ZdMtLIg7lpe7Y$ 
>>>>>>>
>>>>>>> Thanks,
>>>>>>>
>>>>>>> Matt
>>>>>>>
>>>>>>>> "point X has a number of DoF not divisible by 2 field components"
>>>>>>>>
>>>>>>>> Is this a hint that I should have 3 x dim new DoF, and simply not deal with entries that I don't need? Or I am not setting up the section properly (works with just one field)?
>>>>>>>>
>>>>>>>> Thank you,
>>>>>>>> Noam.
>>>>>>>
>>>>>>> --
>>>>>>>
>>>>>>> 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://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/*5D(http:/*www.cse.buffalo.edu/*knepley/)__;fiUvfg!!G_uCfscf7eWS!ZJz_pZR7M_mMi6gIDNMkNQy_6DO0K3sPQyCW-lCG3GWkYxEDauwvsaPpgiykvB00ywn3uvdURZsfZd4cNU2ZdMtLIofCHV1h$ 
>>>>>
>>>>> --
>>>>>
>>>>> 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://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/*5D(http:/*www.cse.buffalo.edu/*knepley/)__;fiUvfg!!G_uCfscf7eWS!ZJz_pZR7M_mMi6gIDNMkNQy_6DO0K3sPQyCW-lCG3GWkYxEDauwvsaPpgiykvB00ywn3uvdURZsfZd4cNU2ZdMtLIofCHV1h$ 
>>>
>>> --
>>>
>>> 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://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/*5D(http:/*www.cse.buffalo.edu/*knepley/)__;fiUvfg!!G_uCfscf7eWS!ZJz_pZR7M_mMi6gIDNMkNQy_6DO0K3sPQyCW-lCG3GWkYxEDauwvsaPpgiykvB00ywn3uvdURZsfZd4cNU2ZdMtLIofCHV1h$ 
>
> --
>
> 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://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/*5D(http:/*www.cse.buffalo.edu/*knepley/)__;fiUvfg!!G_uCfscf7eWS!ZJz_pZR7M_mMi6gIDNMkNQy_6DO0K3sPQyCW-lCG3GWkYxEDauwvsaPpgiykvB00ywn3uvdURZsfZd4cNU2ZdMtLIofCHV1h$ 
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20260508/1f7afab5/attachment-0001.html>