Hi Klaus,
you can do assembling as usual for cells, just using FEFaceFalues, (snippet
from
my code):
FEValues<dim> fe_values (dof_handler.get_fe(), quadrature,
//setup FE values
UpdateFlags(update_values |
update_gradients |
update_q_points |
update_JxW_values));
FEFaceValues<dim> fe_face_values (dof_handler.get_fe(), face_quadrature,
//and FE face values
UpdateFlags(update_values |
update_gradients |
update_q_points |
update_normal_vectors |
update_JxW_values));
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
if ((cell->subdomain_id() == subdomain_id) || (subdomain_id ==
types::invalid_subdomain_id))
{
cell_rhs = 0;
fe_values.reinit (cell);
rhs_functions->body_value_list (fe_values.get_quadrature_points(),
//get the body force value
body_rhs_values);
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const unsigned int component_i =
dof_handler.get_fe().system_to_component_index(i).first;
cell_rhs(i) += (body_rhs_values[q_point](component_i) *
//and integrate, this is (f,v)_cell
fe_values.shape_value(i,q_point) *
fe_values.JxW(q_point));
};
for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell;
++face) //loop over cell faces
if (cell->face(face)->at_boundary()) //if
we are on the neumann boundary
if (cell->face(face)->boundary_indicator() !=
BoundaryIndicators::dirichlet_boundary)
{
fe_face_values.reinit (cell,face);
rhs_functions->boundary_value_list
(fe_face_values.get_quadrature_points(), //get the boundary values
fe_face_values.get_normal_vectors(),
boundary_rhs_values);
for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const unsigned int component_i =
dof_handler.get_fe().system_to_component_index(i).first;
cell_rhs(i) += boundary_rhs_values[q_point](component_i)
* //and integrate, this is (g,v)_face
fe_face_values.shape_value(i, q_point) *
fe_face_values.JxW(q_point);
};
};
cell->get_dof_indices (local_dof_indices);
constraints.distribute_local_to_global (cell_rhs,
local_dof_indices,
rhs);
};
May this helps,
Best,
Martin
________________________________
Von: Klaus Fischer <[email protected]>
An: [email protected]
Gesendet: Dienstag, den 5. April 2011, 14:24:10 Uhr
Betreff: Re: [deal.II] Handle boundary terms
Hi,
I think we're talking about different things here. I'm totally aware of the
boundary conditions of my system and found examples how to implement these
constraints using deal.II. What I meant was looking for example at the
"vector-valued problems" module description, one has the final form (v,u) -
(div
v,p) - (q, div u) = (q, f) + (nv, p)_\Omega.
In the implementation however the last term in the above form, the boundary
term, is neglected. I'm now wondering how would this look like if I don't
neglect it. Somehow I then have to be able to evaluate the shape functions and
their gradients on the surface of each cell. How can this be done?
Best regards,
Klaus
___________________________________________________________
Schon gehört? WEB.DE hat einen genialen Phishing-Filter in die
Toolbar eingebaut! http://produkte.web.de/go/toolbar
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
