Hi,
I need some help here ...
I am having trouble to understand what is going on with a boundary integral.
I have solved a Poisson problem for pressure in a 2D mesh, as a model for a
incompressible water flow in porous media.
The mesh (Tri6 elements) has a near-circular polygon in the middle
representing an injection well.
RHS of the equations are zero except in boundaries dofs. Dirichlet boundary
conditions are imposed using a penalty both on the well and on the external
boundaries to represent fixed pressure. (the solution looks ok - seems to
work)
To validate the solution, I am trying to estimate the water flow both in
the well and in the external boundaries (they must be identical and match
the analytical solution). I therefore calculate the normal pressure
gradient and integrate in the whole boundary path.
When I use first order variable, the results are "almost correct". But when
I use second order , it is absolutely off.
I am using the code below to do the integration. I tried many alternatives,
no success.
Any idea of what I am doing wrong?
Thanks in advance,
Renato
//// Compute boundary integral (internal - well ; and external -limit of
problem domain)
const DofMap & dof_map = _system.get_dof_map();
FEType fe_type = dof_map.variable_type(0);
const unsigned int dim = _lm_mesh.mesh_dimension();
UniquePtr<FEBase> fe_face (FEBase::build(dim, fe_type));
const vector<vector<Real>> & phi_face = fe_face->get_phi();
const vector<vector<RealGradient>> & dphi_face = fe_face->get_dphi();
const vector<Real> & JxW_face = fe_face->get_JxW();
const vector<Point> & normals = fe_face->get_normals();
QGauss qface(dim-1, FIFTH);
fe_face->attach_quadrature_rule (&qface);
const libMesh::BoundaryInfo & boundary_info =
_lm_mesh.get_boundary_info();
MeshBase::const_element_iterator el =
_lm_mesh.active_elements_begin();
const MeshBase::const_element_iterator end_el =
_lm_mesh.active_elements_end();
vector<Number> dpn; dpn.push_back(0); dpn.push_back(0);
int conta0 = 0, conta1 = 0;
for ( ; el != end_el ; ++el)
{
const Elem * elem = *el;
vector<dof_id_type> dof_indices;
dof_map.dof_indices (elem, dof_indices);
const unsigned int num_dofs = cast_int<unsigned int>
(dof_indices.size());
for (unsigned int side=0; side<elem->n_sides(); side++) {
vector<boundary_id_type> vec;
boundary_info.boundary_ids( elem, side, vec );
vector<boundary_id_type>::const_iterator it=vec.begin();
vector<boundary_id_type>::const_iterator it_end=vec.end();
for ( ; it!=it_end; it++ ) {
unsigned int btype = *it; // btype=0:inner(well) ;
btype=1:outer(external limits)
fe_face->reinit(elem, side);
double grad_int_norm = 0;
for (unsigned int qp=0; qp<fe_face->n_quadrature_points(); qp++) {
Gradient grad_at_qn = 0;
for (std::size_t i=0; i<dphi_face.size(); i++)
grad_at_qn += _sys_solution[dof_indices[i]] * dphi_face[i][qp];
Gradient grad_at_qn_jw = JxW_face[qp] * grad_at_qn;
grad_int_norm += grad_at_qn_jw * normals[qp];
}
dpn[btype] += grad_int_norm;
}
}
} // foreach elem
cout << dpn[0]; // well boundary integration result
cout << dpn[1]; // external boundary integration result
------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
_______________________________________________
Libmesh-users mailing list
Libmesh-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-users
