Hello Barbel, thanks for your reply.

What I'm trying to solve is a heat conduction problem in a three
dimensional cube with a hole in the middle (representing a pipe
transporting fluid, a coupled 1d advective problem).

Here is a draft of the 3d domain with the boundary indicators:

        ^z /                            /
| /                            /
|/             1              /
------------------------------
|                            |
|                            |
|                            |
|                            |
| 0        <circular         | 0
|          hole here         |
|          with b.i. 3>      |
|                            |  /
|                            | /
|                            |/
------------------------------ --> x
       /        2
      /
     /y

+Boundary indicator 0 is for heat flux zero. Neumann boundary condition
equal to 0. In the code I left this boundary unspecified because that's
the way to define a neumann boundary value equal to zero, is that right?

+Boundary indicator 1 is for robin boundary condition. The top surface
is under convective heat flux with an external constant temperature of
25 and a convective coefficient of 1. I define this boundary in the
assembly step as:

for (; cell!=endc; ++cell)
.....
  for (unsigned int face = 0; face < GeometryInfo<3>::faces_per_cell;
++face)
  {
   if (cell->face(face)->at_boundary()
       &&
       cell->face(face)->boundary_indicator() == 1)
      {
      fe_face_values.reinit (cell,face);
   for (unsigned int q_face_point = 0; q_face_point < n_face_q_points;
++q_face_point)
        {
         for (unsigned int i = 0; i < dofs_per_cell; ++i)
          {
            for (unsigned int j = 0; j < dofs_per_cell; ++j)
            {
              cell_matrix (i,j) += (1 *
      fe_face_values.shape_value(i,q_face_point) *
      fe_face_values.shape_value(j,q_face_point) *
      fe_face_values.JxW(q_face_point)); 
      }
      cell_rhs (i) += (1 * /*constant temperature*/ 25 *
          fe_face_values.shape_value(i,q_face_point) *
    fe_face_values.JxW(q_face_point));
          }
        }
      }
   }
.....

+Boundary indicator 3 is for robin boundary condition (the hole's
interior boundary). I define this boundary in the assembly step as
(please note how I'm asking for the values from the 1d advective
problem):

//this loop comes immediately after the previous one
for (unsigned int face = 0; face < GeometryInfo<3>::faces_per_cell;
++face)
{
   if (cell->face(face)->at_boundary()
       &&
       cell->face(face)->boundary_indicator() == 3)
      {
        fe_face_values.reinit (cell,face);
        for (unsigned int q_face_point = 0; q_face_point <
n_face_q_points; ++q_face_point)
      {
      for (unsigned int i = 0; i < dofs_per_cell; ++i)
         {
          for (unsigned int j = 0; j < dofs_per_cell; ++j)
           {
              cell_matrix (i,j) += (1 *
fe_face_values.shape_value(i,q_face_point) *
fe_face_values.shape_value(j,q_face_point) *
fe_face_values.JxW(q_face_point));
      }
          cell_rhs (i) += (1 *
   /*solution from 1d problem using y coordinate in 3d problem*/
                VectorTools::point_value (dof_handler_1d,
solution_1d,Point<1> (cell->face(face)->center()(1))) *
fe_face_values.shape_value (i,q_face_point) * 
fe_face_values.JxW(q_face_point)); 
          }
          }
}
       }
....

+ and finally, boundary indicator 2. This is for the bottom face of my
cube. It is defined as a Dirichlet boundary condition (constant
temperature of 10C). I think it is necessary to define at least one
Dirichlet boundary condition for the Laplace equation to be solvable, is
this right? I define it as:

std::map<unsigned int,double> boundary_values;
VectorTools::interpolate_boundary_values (dof_handler_3d,
      2,
      ConstantFunction<3>(10),
      boundary_values);

MatrixTools::apply_boundary_values (boundary_values,
        system_matrix_3d,
        solution_3d,
        system_rhs_3d);




So, my question would be how do I tell the Kelly Error Estimator which
are the boundaries with Neumann boundary conditions and how I specify
its values? (for the current problem I'm not using any function for
Neumann's boundaries and I actually giving them a zero value). And, do I
need to give information about the robin boundary conditions to the
Kelly Error Estimator? If so, how do I do this? I'm not using any
function.

I hope my question is a bit more clear now. Perhaps I have misunderstood
completely how to use the Kelly Error Estimator for generating adaptive
meshes. 



Kind Regards
Javier Munoz
--- Begin Message ---
Dear Javier,

in step-6 only Dirichlet boundary conditions are imposed. So if you would like to use Neumann boundary conditions, you can't you the same functions. For this kind of boundary conditions you would have to implement them in a weak sense (Nitsche's trick).

If you really want Dirichlet conditions and you have two boundary indicators and different boundary functions, you can handle this as

std::map<unsigned int,double> boundary_values;

VectorTools::interpolate_boundary_values (dof_handler,
                                          1,
                                          bundary_function_1(),
                                          boundary_values);

VectorTools::interpolate_boundary_values (dof_handler,
                                          2,
                                          bundary_function_2(),
                                          boundary_values);

MatrixTools::apply_boundary_values (boundary_values,
                                       system_matrix,
                                       solution,
                                       system_rhs);

Does this help? If not, please let us know!

Best,
Bärbel

--
Bärbel Janssen
Mathematisches Institut
Universität Bern

Sidlerstrasse 5, CH-3012 Bern
Tel: +41 (0)31 631 88 47
[email protected]

On Mon, 7 May 2012, Javier Muñoz wrote:

Hello all


I'm trying to use the adaptive refinement capabilities of deal.ii. I was
reading step 6, where adaptivity is introduced for first time.

I have two questions:

The kelly error estimator expects a map with relations between boundary
indicators and corresponding objects describing neumann boundary values.
So, for what I understand if I have, lets say, 2 boundaries with
indicators '1' and '2', I need to put them in the map function with the
respective functions that gave me the values for these boundaries.
Something like:

'1' boundary_function_1()
'2' boundary_function_2()

How do I define this map function?

And, what happen if I use robin boundary conditions?  Do I need to put
them in the map function? If so, how do I do it?

Finally, I would like to expose what is the problem I'm trying to solve
because I'm not sure if the way I'm doing it the more efficient one
(probably not). Any advice would be highly appreciated.

The problem is a heat pipe. I have a 3d cubic domain with a cylindrical
hole. The problem equation in the cube domain is heat diffusion and is
subjected to robin and neumann boundary conditions on the external
surfaces and a robin boundary condition on the hole surface. This
boundary condition depends on the solution of a 1d advection problem.
The 1d domain has an internal heat generation depending on the values on
the internal hole surface in the 3d domain.
I iterate the assembling and solve steps until the difference in both
norm solutions is less than .5.


Kind regards
Javier Munoz




_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii

--- End Message ---
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii

Reply via email to