Hello
I am implementing meshworker into step-33. I have defined class Integrator
to contain the MeshWorker and the integrate functions. When
MeshWorker::integration_loop is called, I am getting segmentation error.
This happens before any of the integrate_* functions are called (I have put
cout statements at beginning of each integrate_* function and these are not
executed). Since the segmentation fault is inside the integration_loop, I am
not able to find the mistake. I have attached some code sections below.
Thanks
praveen
template <int dim>
class Integrator
{
public:
Integrator (const dealii::DoFHandler<dim>& dof_handler)
:
dof_info (dof_handler) {};
dealii::MeshWorker::IntegrationInfoBox<dim> info_box;
dealii::MeshWorker::DoFInfo<dim> dof_info;
dealii::MeshWorker::Assembler::SystemSimple<dealii::TrilinosWrappers::SparseMatrix,
dealii::Vector<double> > assembler;
typedef dealii::MeshWorker::DoFInfo<dim> DoFInfo;
typedef dealii::MeshWorker::IntegrationInfo<dim> CellInfo;
static void integrate_cell_term (const dealii::Vector<double>&
current_solution,
const dealii::Vector<double>&
old_solution,
const Parameters::AllParameters<dim>&
parameters,
DoFInfo& dinfo, CellInfo& info);
static void integrate_boundary_term (const dealii::Vector<double>&
current_solution,
const dealii::Vector<double>&
old_solution,
constParameters::AllParameters<dim>& parameters,
DoFInfo& dinfo, CellInfo& info);
static void integrate_face_term (const dealii::Vector<double>&
current_solution,
const dealii::Vector<double>&
old_solution,
const Parameters::AllParameters<dim>&
parameters,
DoFInfo& dinfo1, DoFInfo& dinfo2,
CellInfo& info1, CellInfo& info2);
};
//------------------------------------------------------------------------------
// Create mesh worker for integration
//------------------------------------------------------------------------------
template <int dim>
void ConservationLaw<dim>::setup_mesh_worker (Integrator<dim>& integrator)
{
std::cout << "Setting up mesh worker ...\n";
MeshWorker::IntegrationInfoBox<dim>& info_box = integrator.info_box;
MeshWorker::DoFInfo<dim>& dof_info = integrator.dof_info;
MeshWorker::Assembler::SystemSimple<TrilinosWrappers::SparseMatrix,
Vector<double> >&
assembler = integrator.assembler;
const unsigned int n_gauss_points = dof_handler.get_fe().degree+1;
info_box.initialize_gauss_quadrature(n_gauss_points,
n_gauss_points,
n_gauss_points);
info_box.initialize_update_flags ();
info_box.add_update_flags_all (update_values |
update_quadrature_points |
update_JxW_values);
info_box.add_update_flags_cell (update_gradients);
info_box.add_update_flags_boundary (update_normal_vectors);
info_box.add_update_flags_face (update_normal_vectors);
info_box.initialize (fe, mapping);
assembler.initialize (system_matrix, right_hand_side);
}
template <int dim>
void ConservationLaw<dim>::assemble_system (Integrator<dim>& integrator)
{
std::cout << "start assemble\n";
MeshWorker::integration_loop<dim, dim>
(dof_handler.begin_active(), dof_handler.end(),
integrator.dof_info,
integrator.info_box,
boost::bind(&Integrator<dim>::integrate_cell_term,
current_solution,
old_solution,
parameters,
_1, _2),
boost::bind(&Integrator<dim>::integrate_boundary_term,
current_solution,
old_solution,
parameters,
_1, _2),
boost::bind(&Integrator<dim>::integrate_face_term,
current_solution,
old_solution,
parameters,
_1, _2, _3, _4),
integrator.assembler, true);
system_matrix.compress ();
}
Then this is the solution part where newton method is used:
unsigned int nonlin_iter = 0;
current_solution = predictor;
Integrator<dim> integrator (dof_handler);
setup_mesh_worker (integrator);
while (true)
{
system_matrix = 0;
right_hand_side = 0;
assemble_system (integrator); etc...
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
