[petsc-users] Additiional DoF per cell

Matthew Knepley knepley at gmail.com
Mon May 11 06:41:57 CDT 2026 

On Mon, May 11, 2026 at 5:25 AM Noam T. <dontbugthedevs at proton.me> wrote:

>
> I have a field that gives me five constant values in every cell
>
> Then you have a P0 field with five components, which will add 5 dofs to
> the DS.
>
>
> Alright, thanks for the clarification.
>
> One additional question: for computing the jacobian/residual, I'm using  DMSNESSetFunction/Jacobian(Local).
> It seems using this function expects a jacobian of full size (for all
> fields), so e.g. 11x11 (8 for u field, 3 for p field). Is that correct? At
> least when looking at MatView after AssemblyBegin/End, that's what I get
> (minus the number of constrained DoF, if any).
> Similarly for the residual.
>

Yes, that is completely true.

More explanation. For the residual, you can always use DMCreateSubDM() to
compute the residual over a subset of fields. However, Plex always computes
square Jacobians. I did write the code to compute rectangular pieces of the
Jacobian, but no one ever used it and eventually removed for maintenance
reasons.

  Thanks,

     Matt


> Thanks,
> Noam
>
>
> On Friday, May 8th, 2026 at 10:29 PM, Matthew Knepley <knepley at gmail.com>
> wrote:
>
> On Fri, May 8, 2026 at 3:34 PM Noam T. <dontbugthedevs at proton.me> wrote:
>
>> I was indeed confused before regarding how the element ended up with 3
>> unknowns.
>>
>> So for a P1, having three unknowns, means the field is interpolated as p
>> = po + p1*x + p2*y?
>> Then P2 has 6 components (total dim = 14), which agrees with a complete
>> 2nd order polynomial.
>>
>> Then my question is, how can I have any other number of components? In
>> particular, some of the implementations I am looking at use 4 and 7
>> additional DoF. Will this require a new field per DoF?
>>
>
> No, you start with how these fields are supposed to behave. If you say
>
> I have a field that gives me five constant values in every cell
>
> Then you have a P0 field with five components, which will add 5 dofs to
> the DS. If instead you say
>
> I have a field that gives two values linearly interpolated over each cell
>
> Then you have a P1 field with two components, which will add 6 dofs (with
> triangles) to the DS.
>
> Thanks,
>
> Matt
>
>> Thanks,
>> Noam
>> On Friday, May 8th, 2026 at 7:13 PM, Matthew Knepley <knepley at gmail.com>
>> wrote:
>>
>> On Fri, May 8, 2026 at 2:54 PM Noam T. <dontbugthedevs at proton.me> wrote:
>>
>>> Code attached.
>>>
>>
>> I think you misunderstand the definition of the element. I do not think
>> you want 3 components for your pressure. You want 1 component. It is a
>> scalar variable. In my original example, I had a single component.
>> Components create tensor spaces.
>>
>> If you are using discontinuous P1, then it has 3 unknowns, which is what
>> I think you are trying to specify. You should change the space you are
>> using if you want to change the number of variables.
>>
>> Does that clear it up?
>>
>> Thanks,
>>
>> Matt
>>
>>> Noticed from PetscSectionView:
>>>
>>> PetscSectionSym Object: 1 MPI process
>>> type: label
>>> Label 'depth'
>>> Symmetry for stratum value 0 (0 dofs per point): no symmetries
>>> Symmetry for stratum value 1 (0 dofs per point): no symmetries
>>> Symmetry for stratum value 2 (12 dofs per point):
>>> <---------------------- 12 DoF ?
>>> Orientation range: [-4, 4)
>>> Symmetry for stratum value -1 (0 dofs per point): no symmetries
>>>
>>> Thanks,
>>> Noam
>>>
>>> On Friday, May 8th, 2026 at 3:58 PM, Noam T. <dontbugthedevs at proton.me>
>>> wrote:
>>>
>>> This is likely to be a Section construction problem again. The Section
>>> fields are separated by DMPlexVecGetClosure(), sp you must not have defined
>>> the fields in the Section. If you put these calls into the code that I sent
>>> you, what do you get?
>>>
>>>
>>> If I use the exact same arrays/values for constructing the section, I
>>> get the same closure arrangement in the code you provided. Not surprising I
>>> guess? A seemingly random number of coordinates / zeroes / more coordinates.
>>>
>>> I checked example 14 (
>>> https://urldefense.us/v3/__https://petsc.org/release/src/dm/impls/plex/tutorials/ex14.c.html__;!!G_uCfscf7eWS!blgAbufRlZ05shY0R_LZ-jc4h4ItGucyIJbrAy_Gsds7p62MI2Fpg_5uf4wgmC6fJwAr0yCf8d3KucejXyYH$ ) for
>>> constructing the section. These are the values I used:
>>>
>>> - n_fields = 2
>>>
>>> - numDof [n_fields*(dim+1)] = numDof[2*(2+1)] = [2 0 0 0 0 N]
>>>
>>> The u field is defined in the nodes, with 2 components (u, v), hence the
>>> first "2". Then no components for edges/faces.
>>> For the p field, what I am trying to simulate is having nc (independent)
>>> DoF, sort of like having nc scalar fields, if that makes sense? Hence using:
>>>
>>> - N = nc
>>> - Also using the flag "-pres_petscspace_degree 0"
>>>
>>> But I'm not sure about those values.
>>>
>>> Then for simplicity numBC = 0 and all other arguments are NULL (will
>>> deal with boundary conditions later).
>>>
>>> With these arguments, I get the same as before e.g.
>>>
>>> -- nc = 3 (total dimension = 11)
>>> closure = y1 x2 y2 x3 y3 0 0 0 x0 y0 x1
>>>
>>> Thanks,
>>> Noam
>>>
>>> On Friday, May 8th, 2026 at 1:25 PM, Matthew Knepley <knepley at gmail.com>
>>> wrote:
>>>
>>> On Fri, May 8, 2026 at 8:01 AM Noam T. <dontbugthedevs at proton.me> wrote:
>>>
>>>> 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!blgAbufRlZ05shY0R_LZ-jc4h4ItGucyIJbrAy_Gsds7p62MI2Fpg_5uf4wgmC6fJwAr0yCf8d3KuZIyGXkb$ ), 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.
>>>>
>>>
>>> This is likely to be a Section construction problem again. The Section
>>> fields are separated by DMPlexVecGetClosure(), sp you must not have defined
>>> the fields in the Section. If you put these calls into the code that I sent
>>> you, what do you get?
>>>
>>> Thanks,
>>>
>>> Matt
>>>
>>>> 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__;!!G_uCfscf7eWS!blgAbufRlZ05shY0R_LZ-jc4h4ItGucyIJbrAy_Gsds7p62MI2Fpg_5uf4wgmC6fJwAr0yCf8d3KuZIyGXkb$ 
>>>>>>>> <https://urldefense.us/v3/__https://petsc.org/release/src/snes/tutorials/ex77.c.html__;!!G_uCfscf7eWS!aZNOxZZeoHYqe6lRZV0wHMVVTu3YIEfc1Sr-dca7xiGVfG3enULFD4_g6_fxWIa7A_Cu9LcXU-MuxRj4sOwDRyWqvsdGdrEd$>)
>>>>>>>> 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!blgAbufRlZ05shY0R_LZ-jc4h4ItGucyIJbrAy_Gsds7p62MI2Fpg_5uf4wgmC6fJwAr0yCf8d3KuSgpU3xx$ 
>>>>>>>
>>>>>>> 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/__;fg!!G_uCfscf7eWS!blgAbufRlZ05shY0R_LZ-jc4h4ItGucyIJbrAy_Gsds7p62MI2Fpg_5uf4wgmC6fJwAr0yCf8d3KuTKxFdNr$ 
>>>>>>> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!blgAbufRlZ05shY0R_LZ-jc4h4ItGucyIJbrAy_Gsds7p62MI2Fpg_5uf4wgmC6fJwAr0yCf8d3KucmJd5wz$ >
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>
>>>>>> --
>>>>>> 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/__;fg!!G_uCfscf7eWS!blgAbufRlZ05shY0R_LZ-jc4h4ItGucyIJbrAy_Gsds7p62MI2Fpg_5uf4wgmC6fJwAr0yCf8d3KuTKxFdNr$ 
>>>>>> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!blgAbufRlZ05shY0R_LZ-jc4h4ItGucyIJbrAy_Gsds7p62MI2Fpg_5uf4wgmC6fJwAr0yCf8d3KucmJd5wz$ >
>>>>>>
>>>>>>
>>>>>>
>>>>>
>>>>> --
>>>>> 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/__;fg!!G_uCfscf7eWS!blgAbufRlZ05shY0R_LZ-jc4h4ItGucyIJbrAy_Gsds7p62MI2Fpg_5uf4wgmC6fJwAr0yCf8d3KuTKxFdNr$ 
>>>>> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!blgAbufRlZ05shY0R_LZ-jc4h4ItGucyIJbrAy_Gsds7p62MI2Fpg_5uf4wgmC6fJwAr0yCf8d3KucmJd5wz$ >
>>>>>
>>>>>
>>>>>
>>>>
>>>> --
>>>> 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/__;fg!!G_uCfscf7eWS!blgAbufRlZ05shY0R_LZ-jc4h4ItGucyIJbrAy_Gsds7p62MI2Fpg_5uf4wgmC6fJwAr0yCf8d3KuTKxFdNr$ 
>>>> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!blgAbufRlZ05shY0R_LZ-jc4h4ItGucyIJbrAy_Gsds7p62MI2Fpg_5uf4wgmC6fJwAr0yCf8d3KucmJd5wz$ >
>>>>
>>>>
>>>>
>>>
>>> --
>>> 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/__;fg!!G_uCfscf7eWS!blgAbufRlZ05shY0R_LZ-jc4h4ItGucyIJbrAy_Gsds7p62MI2Fpg_5uf4wgmC6fJwAr0yCf8d3KuTKxFdNr$ 
>>> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!blgAbufRlZ05shY0R_LZ-jc4h4ItGucyIJbrAy_Gsds7p62MI2Fpg_5uf4wgmC6fJwAr0yCf8d3KucmJd5wz$ >
>>>
>>>
>>>
>>>
>>
>> --
>> 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/__;fg!!G_uCfscf7eWS!blgAbufRlZ05shY0R_LZ-jc4h4ItGucyIJbrAy_Gsds7p62MI2Fpg_5uf4wgmC6fJwAr0yCf8d3KuTKxFdNr$ 
>> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!blgAbufRlZ05shY0R_LZ-jc4h4ItGucyIJbrAy_Gsds7p62MI2Fpg_5uf4wgmC6fJwAr0yCf8d3KucmJd5wz$ >
>>
>>
>>
>
> --
> 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/__;fg!!G_uCfscf7eWS!blgAbufRlZ05shY0R_LZ-jc4h4ItGucyIJbrAy_Gsds7p62MI2Fpg_5uf4wgmC6fJwAr0yCf8d3KuTKxFdNr$ 
> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!blgAbufRlZ05shY0R_LZ-jc4h4ItGucyIJbrAy_Gsds7p62MI2Fpg_5uf4wgmC6fJwAr0yCf8d3KucmJd5wz$ >
>
>
>

-- 
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/__;fg!!G_uCfscf7eWS!blgAbufRlZ05shY0R_LZ-jc4h4ItGucyIJbrAy_Gsds7p62MI2Fpg_5uf4wgmC6fJwAr0yCf8d3KuTKxFdNr$  <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!blgAbufRlZ05shY0R_LZ-jc4h4ItGucyIJbrAy_Gsds7p62MI2Fpg_5uf4wgmC6fJwAr0yCf8d3KucmJd5wz$ >
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20260511/2f295ce4/attachment-0001.html>

More information about the petsc-users mailing list