> I forgot that earlier. Could you come up with a small program that shows
> the
> segfault? We should have caught the wrong vector size via an assertion,
> and
> I'd like to add that assertion.
>
Certainly--please find attached.
To clarify, deal.II *did* throw some type of exception, I just wasn't able
to interpret it and opted to use GDB instead; I'm a Python programmer and
not a C++ expert by any stretch of the imagination :)
The program terminates with:
An error occurred in line <3196> of file
</home/foucartc/software/dealii-9.2.0/source/fe/fe_values.cc> in function
void dealii::internal::do_function_values(const typename
VectorType::value_type*, const dealii::Table<2, double>&, const
dealii::FiniteElement<dim, spacedim>&, const std::vector<unsigned int>&,
dealii::ArrayView<Number>, bool, unsigned int) [with int dim = 2; int
spacedim = 2; VectorType = dealii::Vector<double>; typename
VectorType::value_type = double]
The violated condition was:
static_cast<typename ::dealii::internal::argument_type<void( typename
std::common_type<decltype(values.size()),
decltype(n_quadrature_points)>::type)>::type>(values.size()) ==
static_cast<typename ::dealii::internal::argument_type<void( typename
std::common_type<decltype(values.size()),
decltype(n_quadrature_points)>::type)>::type>(n_quadrature_points)
Additional information:
Dimension 0 not equal to 3.
Stacktrace:
-----------
#0
/home/foucartc/CorbinFoucart_repos/dealii_applications/lib/libdeal_II.g.so.9.2.0:
void dealii::internal::do_function_values<2, 2, dealii::Vector<double>
>(dealii::Vector<double>::value_type const*, dealii::Table<2, double>
const&, dealii::FiniteElement<2, 2> const&, std::vector<unsigned int,
std::allocator<unsigned int> > const&,
dealii::ArrayView<dealii::Vector<double>, dealii::MemorySpace::Host>, bool,
unsigned int)
#1
/home/foucartc/CorbinFoucart_repos/dealii_applications/lib/libdeal_II.g.so.9.2.0:
void dealii::FEValuesBase<2, 2>::get_function_values<dealii::Vector<double>
>(dealii::Vector<double> const&,
std::vector<dealii::Vector<dealii::Vector<double>::value_type>,
std::allocator<dealii::Vector<dealii::Vector<double>::value_type> > >&)
const
#2 ./step: DGVelocityField<2>::access_interface_values()
#3 ./step: main
--------------------------------------------------------
Aborted (core dumped)
As for implementing FEInterfaceValues::get_function_jumps et al: I have a
deadline by which I'd really like these set of codes to be working, so I
can't work on this in the next few weeks. However, I will have to implement
that functionality anyway, so I wouldn't mind working on this sometime in
July.
Best,
Corbin
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/1bfe5faa-2f0e-4b8a-a0b6-015c73aaa469n%40googlegroups.com.
// unification to a depth / surface mesh class
//
#include <iostream>
#include <fstream>
// deal ii imports
#include <deal.II/base/function_lib.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_interface_values.h
using namespace dealii;
template<int dim>
struct DGData
{
FEValues<dim> fe_values;
FEFaceValues<dim> fe_face_values;
FEInterfaceValues<dim> fe_interface_values;
// stores mapping between faces and the cell interior basis
// functions which have support on each face
std::vector<std::vector<unsigned int>> fe_local_support_on_face;
DGData(
const FiniteElement<dim> &fe,
const Quadrature<dim> &cell_quadrature,
const Quadrature<dim-1> &face_quadrature,
const UpdateFlags cell_update_flags = update_values
| update_gradients
| update_quadrature_points
| update_JxW_values,
const UpdateFlags interface_update_flags = update_values
| update_quadrature_points
| update_JxW_values
| update_gradients
| update_normal_vectors)
: fe_values(fe, cell_quadrature, cell_update_flags)
, fe_face_values(fe, face_quadrature, interface_update_flags)
, fe_interface_values(fe, face_quadrature, interface_update_flags)
, fe_local_support_on_face(GeometryInfo<dim>::faces_per_cell)
{ }
};
template <int dim>
class DGVelocityField
{
public:
DGVelocityField(
const unsigned int degree,
Triangulation<dim> &grid)
: triangulation(grid)
, fe_sys(FE_DGQ<dim>(degree), dim)
, dofh(triangulation)
, quadrature(degree + 1)
, face_quadrature(degree + 1)
, dgd(fe_sys, quadrature, face_quadrature)
{
// system DOF
dofh.distribute_dofs(fe_sys);
velocity_field.reinit(dofh.n_dofs());
}
Triangulation<dim> &triangulation;
FESystem<dim> fe_sys;
DoFHandler<dim> dofh;
Vector<double> velocity_field;
const QGauss<dim> quadrature;
const QGauss<dim - 1> face_quadrature;
DGData<dim> dgd;
void access_interface_values();
};
template <int dim>
void DGVelocityField<dim>::access_interface_values()
{
for (const auto &cell : dofh.active_cell_iterators())
{
dgd.fe_values.reinit(cell);
for (const auto face_index : GeometryInfo<dim>::face_indices())
{
if (!cell->face(face_index)->at_boundary())
{
dgd.fe_interface_values.reinit(cell,
face_index,
numbers::invalid_unsigned_int,
cell->neighbor(face_index),
cell->neighbor_of_neighbor(face_index),
numbers::invalid_unsigned_int);
// attempt to access neighboring cell's velocities at qpts, SEGFAULT
// assignment should instead involve:
// const unsigned int n_face_qpts = dgd.fe_interface_values.n_quadrature_points;
// std::vector<Vector<double>> vel_vals(n_face_qpts, Vector<double>(dim));
std::vector<Vector<double>> vel_vals;
dgd.fe_interface_values.get_fe_face_values(1).get_function_values(velocity_field, vel_vals);
}
}
} // end cell loop
}
int main()
{
const unsigned int dim = 2;
Triangulation<dim> triangulation;
GridGenerator::subdivided_hyper_cube(triangulation, 2, -1, 1); // 2^dim elements
DGVelocityField<dim> v(2, triangulation);
v.access_interface_values();
return 0;
}
