Impose values inside a material/ not on the boundary

[Markus Mehnert]

Hi everyone,

Instead of prescribing classical Dirichlet boundary conditions (i.e. the values for the solution field on the boundary) I am trying to impose the values for a specific degree of freedom inside my material. My idea was to iterate over all cells and their faces. Then I use the location of the center of the face to identify the faces, that have the right position. When I loop over all degrees of freedom of the face, I can find the DOFs that correspond to the correct base element, put these into a "boundary_values" vector and then use MatrixTools::apply_boundary_values to apply these values (see the following code snippet).

 typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active();
for (; cell!=dof_handler.end(); ++cell)
{
        for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)

                const Point<dim> face_center = cell->face(f)->center();
        if (face_center[2] == 5.0)
        {
                cell->face(f)->get_dof_indices(local_dof_indices_face,0);
                for(unsigned int i = 0; i < dofs_per_face; ++i)
                {
                        const unsigned int i_group     = fe_face.system_to_base_index(i).first.first;
                        if (i_group == V_dof)
                        {
                                boundary_values[local_dof_indices_face[i]]=0;
                        }
                }

        }
}


MatrixTools::apply_boundary_values (boundary_values,
                system_matrix,
                locally_relevant_solution_update,
                system_rhs,false);

Here I try to apply the value for the electric potential (a scalar value, associated with the V_dof) on the middle plane of a cube with the sidelength 10 (i.e. face_center[2] == 5.0)

This works fine, if the polynomial order of the FE_Q element is 1, however, if I increase the order it seems like the value is not applied to the DOFs between the vertices of the cell (see attached picture).  When I apply the boundary conditions on the surfaces using VectorTools::interpolate_boundary_values I do not have this issue.
Can anyone give me a tip, what might be the problem?

Thank you,
Markus

[Wolfgang Bangerth]

Markus,
how do you define dofs_per_face? Do you adjust it for the polynomial degree?
Best
W.


On 2/11/22 10:51, 'Markus Mehnert' via deal.II User Group wrote:
> *** Caution: EXTERNAL Sender ***

>
> Hi everyone,
>
> Instead of prescribing classical Dirichlet boundary conditions (i.e. the
> values for the solution field on the boundary) I am trying to impose the
> values for a specific degree of freedom inside my material. My idea was to
> iterate over all cells and their faces. Then I use the location of the center
> of the face to identify the faces, that have the right position. When I loop
> over all degrees of freedom of the face, I can find the DOFs that correspond
> to the correct base element, put these into a "boundary_values" vector and
> then use MatrixTools::apply_boundary_values to apply these values (see the
> following code snippet).
>

> / typename hp::DoFHandler<dim>::active_cell_iterator

> cell = dof_handler.begin_active();
> for (; cell!=dof_handler.end(); ++cell)
> {
>         for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
>
>                 const Point<dim> face_center = cell->face(f)->center();
>         if (face_center[2] == 5.0)
>         {
>                 cell->face(f)->get_dof_indices(local_dof_indices_face,0);
>                 for(unsigned int i = 0; i < dofs_per_face; ++i)
>                 {
>                         const unsigned int i_group     =
> fe_face.system_to_base_index(i).first.first;
>                         if (i_group == V_dof)
>                         {
>                                 boundary_values[local_dof_indices_face[i]]=0;
>                         }
>                 }
>
>         }
> }
>
>
> MatrixTools::apply_boundary_values (boundary_values,
>                 system_matrix,
>                 locally_relevant_solution_update,

>                 system_rhs,false);/

>
>
> Here I try to apply the value for the electric potential (a scalar value,
> associated with the V_dof) on the middle plane of a cube with the sidelength

> 10 (i.e. /face_center[2] == 5.0/)

>
> This works fine, if the polynomial order of the FE_Q element is 1, however, if
> I increase the order it seems like the value is not applied to the DOFs
> between the vertices of the cell (see attached picture).  When I apply the
> boundary conditions on the surfaces using

> /VectorTools::interpolate_boundary_values /I do not have this issue.

> Can anyone give me a tip, what might be the problem?
>
> Thank you,
> Markus
>

> --
> The deal.II project is located at http://www.dealii.org/
> <https://nam10.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.dealii.org%2F&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cc4862c0f76fd4db4f40c08d9ed8727c3%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637801987794503970%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=a4ySeNYTtAc94%2BX8v3pc9l99EeGx0NF%2B9vVE4UoUjbM%3D&reserved=0>
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fforum%2Fdealii%3Fhl%3Den&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cc4862c0f76fd4db4f40c08d9ed8727c3%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637801987794503970%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=P1P8ynFMK9a%2BeZk8D5tfng8lAfCb%2FkUpG1jKJOsbiVY%3D&reserved=0>
> ---
> You received this message because you are subscribed to the Google Groups
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to dealii+un...@googlegroups.com
> <mailto:dealii+un...@googlegroups.com>.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/5b935808-c22f-4676-9418-6d7ee30228b9n%40googlegroups.com
> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fmsgid%2Fdealii%2F5b935808-c22f-4676-9418-6d7ee30228b9n%2540googlegroups.com%3Futm_medium%3Demail%26utm_source%3Dfooter&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cc4862c0f76fd4db4f40c08d9ed8727c3%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637801987794503970%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=MMXHnlG4USb1BqZxqfGnG2FZ9WGrLsdIkLeP%2B%2BDfvgo%3D&reserved=0>.


--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Markus Mehnert]

Dear Wolfgang,

Thank you for your quick response. The dofs per face are defined as

const unsigned int dofs_per_face = fe_cell.dofs_per_face;

where fe_cell is the FESystem that consists of

fe_cell(FE_Q<dim> (poly_order), dim,  //Displacement
             FE_Q<dim> (poly_order), 1,          //Electric Potential
             FE_Q<dim> (poly_order+1), dim, //Vectorial Electronic order parameter
             FE_DGPMonomial<dim>(poly_order - 1), 1,  // Dilatation
             FE_DGPMonomial<dim>(poly_order - 1), 1),  // Pressure

 I checked the number of dofs and it seems to be correct.

Best,

Markus

[Wolfgang Bangerth]

On 2/11/22 12:12, 'Markus Mehnert' via deal.II User Group wrote:
> *** Caution: EXTERNAL Sender ***
>

> Dear Wolfgang,
>
> Thank you for your quick response. The dofs per face are defined as /
> /
> /
> /
> /const unsigned int dofs_per_face = fe_cell.dofs_per_face;/
> /
> /
> where fe_cell is the/FESystem/ that consists of

>
> fe_cell(FE_Q<dim> (poly_order), dim,  //Displacement
>              FE_Q<dim> (poly_order), 1,          //Electric Potential
>              FE_Q<dim> (poly_order+1), dim, //Vectorial Electronic order
> parameter
>              FE_DGPMonomial<dim>(poly_order - 1), 1,  // Dilatation
>              FE_DGPMonomial<dim>(poly_order - 1), 1),  // Pressure
>
>  I checked the number of dofs and it seems to be correct.

Then it's time to count how many constraints you expect (on a piece of paper),
how many you actually enter into the constraint object, and how many the
constraints object actually holds :-)

Best
W.

[Markus Mehnert]

Dear Wolfgang,

Thank you for your response. I identified the correct dofs by trial and error which worked (9 dofs for a scalar quantity with quadratic polyorder).  However, when I use "fe_face.system_to_base_index().first.first" to identify, which base element these dofs belong to, it does not at all fit the theory as they should all belong to base 1, i.e. the electric potential (I copied the output from my terminal)

The base group of dof 3 is: 1
The base group of dof 10 is: 1
The base group of dof 17 is: 1
The base group of dof 24 is: 1
The base group of dof 31 is: 1
The base group of dof 41 is: 2
The base group of dof 51 is: 0
The base group of dof 61 is: 2
The base group of dof 71 is: 2

This explains, why my former code did not work, however it raises the question, why the dofs do not belong to the (correct) base element. Do you have any idea?

Thank you and kind regards,

Markus

[Wolfgang Bangerth]

On 2/15/22 07:32, 'Markus Mehnert' via deal.II User Group wrote:
>
> Thank you for your response. I identified the correct dofs by trial and error
> which worked (9 dofs for a scalar quantity with quadratic polyorder).

> However, when I use /"fe_face.system_to_base_index().first.first" /to

> identify, which base element these dofs belong to, it does not at all fit the
> theory as they should all belong to base 1, i.e. the electric potential (I
> copied the output from my terminal)
>
> The base group of dof 3 is: 1
> The base group of dof 10 is: 1
> The base group of dof 17 is: 1
> The base group of dof 24 is: 1
> The base group of dof 31 is: 1
> The base group of dof 41 is: 2
> The base group of dof 51 is: 0
> The base group of dof 61 is: 2
> The base group of dof 71 is: 2
>
> This explains, why my former code did not work, however it raises the
> question, why the dofs do not belong to the (correct) base element. Do you
> have any idea?

Markus -- I have to admit that I'm not sure I can follow the details. Can you
create a small program that only creates the element and outputs what you are
interested in, along with an explanation of what you *think* it should output?