This is my assemble function. I can't see that I'm doing anything wrong but
then again, I'm new to the library.
void assemble_poisson(EquationSystems& es,
const std::string& system_name)
{
libmesh_assert (system_name == "Poisson");
const MeshBase& mesh = es.get_mesh();
const unsigned int dim = mesh.mesh_dimension();
LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>
("Poisson");
const DofMap& dof_map = system.get_dof_map();
FEType fe_type = dof_map.variable_type(0);
AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
QGauss qrule (dim, FIFTH);
fe->attach_quadrature_rule (&qrule);
AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
QGauss qface(dim-1, FIFTH);
fe_face->attach_quadrature_rule (&qface);
const std::vector<Real>& JxW = fe->get_JxW();
const std::vector<Point>& q_point = fe->get_xyz();
const std::vector<std::vector<Real> >& phi = fe->get_phi();
const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
DenseMatrix<Number> Ke;
DenseVector<Number> Fe;
std::vector<unsigned int> dof_indices;
MeshBase::const_element_iterator el =
mesh.active_local_elements_begin();
const MeshBase::const_element_iterator end_el =
mesh.active_local_elements_end();
for ( ; el != end_el ; ++el)
{
const Elem* elem = *el;
dof_map.dof_indices (elem, dof_indices);
fe->reinit (elem);
Ke.resize (dof_indices.size(),
dof_indices.size());
Fe.resize (dof_indices.size());
for (unsigned int qp=0; qp<qrule.n_points(); qp++)
{
for (unsigned int i=0; i<phi.size(); i++)
for (unsigned int j=0; j<phi.size(); j++)
Ke(i,j) +=
JxW[qp]*(q_point[qp](0)*dphi[i][qp]*dphi[j][qp]);
{
const Real x = q_point[qp](0);
//const Real y = q_point[qp](1);
//const Real eps = 1.e-3;
const Real fxy = -2/(q_point[qp](0) * q_point[qp](0));
for (unsigned int i=0; i<phi.size(); i++)
Fe(i) += JxW[qp]*fxy*phi[i][qp];
}
}
{
for (unsigned int side=0; side<elem->n_sides(); side++){
short int bc_id = mesh.boundary_info->boundary_id
(elem,side);
if (elem->neighbor(side) == NULL)
{
const std::vector<std::vector<Real> >& phi_face =
fe_face->get_phi();
const std::vector<Real>& JxW_face = fe_face->get_JxW();
const std::vector<Point >& qface_point =
fe_face->get_xyz();
fe_face->reinit(elem, side);
for (unsigned int qp=0; qp<qface.n_points(); qp++)
{
const Real xf = qface_point[qp](0);
const Real yf = qface_point[qp](1);
const Real penalty = 1.e10;
if (bc_id == 0) { //Dirichlet BC
// The boundary value.
const Real value = 2.0;//exact_solution(xf, yf);
for (unsigned int i=0; i<phi_face.size(); i++)
for (unsigned int j=0; j<phi_face.size();
j++)
Ke(i,j) +=
JxW_face[qp]*penalty*phi_face[i][qp]*phi_face[j][qp];
for (unsigned int i=0; i<phi_face.size(); i++)
Fe(i) +=
JxW_face[qp]*penalty*value*phi_face[i][qp];
}
else if (bc_id == 1) //Neumann BC
{
const Real value = -0.5;//neumann_value(xf, yf,
zf);
for (unsigned int i=0; i<phi.size(); i++)
Fe(i) += JxW[qp]*phi[i][qp]*value;
}
}
}
}
}
dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
system.matrix->add_matrix (Ke, dof_indices);
system.rhs->add_vector (Fe, dof_indices);
}
}
Regards
Ted
------------------------------------------------------------------------------
Let Crystal Reports handle the reporting - Free Crystal Reports 2008 30-Day
trial. Simplify your report design, integration and deployment - and focus on
what you do best, core application coding. Discover what's new with
Crystal Reports now. http://p.sf.net/sfu/bobj-july
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
