[petsc-users] Additiional DoF per cell

Noam T. dontbugthedevs at proton.me
Fri May 8 14:34:03 CDT 2026 

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?

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!aFjajGFLkDiyYC1UgjN8b9QCym74RzcADVVBaCHOtEnbFz0QyrOJWDGfPaj-vfN47fuweu5C_zilFW5N6FsfqFdv2WmI-lhp$ ) 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!aFjajGFLkDiyYC1UgjN8b9QCym74RzcADVVBaCHOtEnbFz0QyrOJWDGfPaj-vfN47fuweu5C_zilFW5N6FsfqFdv2SV3HT11$ ), 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*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!aFjajGFLkDiyYC1UgjN8b9QCym74RzcADVVBaCHOtEnbFz0QyrOJWDGfPaj-vfN47fuweu5C_zilFW5N6FsfqFdv2UyvJen8$ ) 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!aFjajGFLkDiyYC1UgjN8b9QCym74RzcADVVBaCHOtEnbFz0QyrOJWDGfPaj-vfN47fuweu5C_zilFW5N6FsfqFdv2cwxv_v8$ 
>>>>>>>>>>>>
>>>>>>>>>>>> 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!aFjajGFLkDiyYC1UgjN8b9QCym74RzcADVVBaCHOtEnbFz0QyrOJWDGfPaj-vfN47fuweu5C_zilFW5N6FsfqFdv2Q8o7bg-$ 
>>>>>>>>>>
>>>>>>>>>> --
>>>>>>>>>>
>>>>>>>>>> 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!aFjajGFLkDiyYC1UgjN8b9QCym74RzcADVVBaCHOtEnbFz0QyrOJWDGfPaj-vfN47fuweu5C_zilFW5N6FsfqFdv2Q8o7bg-$ 
>>>>>>>>
>>>>>>>> --
>>>>>>>>
>>>>>>>> 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!aFjajGFLkDiyYC1UgjN8b9QCym74RzcADVVBaCHOtEnbFz0QyrOJWDGfPaj-vfN47fuweu5C_zilFW5N6FsfqFdv2Q8o7bg-$ 
>>>>>>
>>>>>> --
>>>>>>
>>>>>> 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!aFjajGFLkDiyYC1UgjN8b9QCym74RzcADVVBaCHOtEnbFz0QyrOJWDGfPaj-vfN47fuweu5C_zilFW5N6FsfqFdv2Q8o7bg-$ 
>>>>
>>>> --
>>>>
>>>> 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!aFjajGFLkDiyYC1UgjN8b9QCym74RzcADVVBaCHOtEnbFz0QyrOJWDGfPaj-vfN47fuweu5C_zilFW5N6FsfqFdv2Q8o7bg-$ 
>
> --
>
> 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!aFjajGFLkDiyYC1UgjN8b9QCym74RzcADVVBaCHOtEnbFz0QyrOJWDGfPaj-vfN47fuweu5C_zilFW5N6FsfqFdv2Q8o7bg-$ 
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20260508/d4b12b7b/attachment-0001.html>

More information about the petsc-users mailing list