Hello, I want to use parts of deal.ii to work with a problem which does not impose any constraints on the boundary. However, I run into some problems with boundary dofs and do not see a way to properly solve them.
I am using deali.ii 6.3.1. The attached code serves to illustrate my problem. It was my first attempt to calculate the discretised representation of a 2D function (the constant function in this case), and output it. The result however is not what I want. It is not constant and at the boundaries the value of the reconstructed function is too small. I believe this is the case, because there are contributions to the value of the boundary dofs are implicitly zero, as the cells for these contributions are not existing. Or posed differently, it is a problem of normalisation, as the boundary dof shape functions have a smaller support (only one quarter for a dof in a corner, one half for a dof on a edge of the domain). I guess that this is usually not a problem, when one uses Dirichlet or Neumann boundary conditions, where the values for these dofs are fixed. My idea of a solution would be to scale the contributions to boundary cells by a factor to compensate for the wrong normalisation. But I do not see how to determine this factor reliably and for all dimensions. I have to find out for all the boundary dofs (which I can obtain with DoFTools::extract_dofs_with_support_on_boundary) how many contributing cells are missing. What is an appropriate way to do that? Thank you in advance Regards Johannes Reinhardt
#include <fstream>
#include <iostream>
#include <grid/tria.h>
#include <grid/tria_accessor.h>
#include <grid/tria_iterator.h>
#include <grid/grid_out.h>
#include <numerics/data_out.h>
#include <fe/fe_q.h>
#include <fe/fe_values.h>
#include <grid/grid_generator.h>
int main(){
//setup grid and fevalues
FE_Q<2> fe(1);
QGauss<2> quad(1);
Triangulation<2> triangulation;
DoFHandler<2> dof_handler(triangulation);
GridGenerator::hyper_cube (triangulation,-1,1);
triangulation.refine_global(2);
dof_handler.distribute_dofs(fe);
FEValues<2> fe_values(fe, quad,update_values | update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quad.size();
std::vector<unsigned int> local_dof_indices (dofs_per_cell);
//the vector containing the discretised function
Vector<double> discrete (dofs_per_cell);
discrete = 0.0;
//loop over cells and sum up the contributions to the coeffs for the dofs
//for the constant function with value 1.0
DoFHandler<2>::active_cell_iterator cell = dofh.begin_active();
DoFHandler<2>::active_cell_iterator endc = dofh.end();
for (; cell!=endc; ++cell){
fe_values.reinit(cell);
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
dst(local_dof_indices[i]) += fe_values.shape_value(i,q_point) *
1.0 * fe_values.JxW(q_point);
}
//output the solution
DataOut<2> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(src,"Constant");
data_out.build_patches();
std::ofstream conv_out("constant.vtk");
data_out.write_vtk(conv_out);
}
<<attachment: notconstant.png>>
_______________________________________________ dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
