> Dear Martin
>
> We do not offer convenience access for all data structures on faces yet,
> so feel free to suggest what you need.
>
1) functions similar to get_cell_iterator for getting iterator to
interior/exterior cell + face index:
std::pair<DoFHandler< dim >::cell_iterator, unstigned int>
get_face_iterator(const unsigned int face_batch_number,
const unsigned int vector_number, const bool iterior=true,
const unsigned int fe_component=0)
2) getter for cell_level_index so that the data from get_face_info could be
used easily (just in case)
> That said, for the access of rho and mu on an element, we provide the
> following access function:
>
>
> https://www.dealii.org/current/doxygen/deal.II/classFEEvaluationBase.html#a442f1ee55e6e80a58e19d53c7184c1a7
>
> This FEEvaluationBase::read_cell_data function expects an
> AlignedVector<VectorizedArray<Number>> as the argument with as many
> components as there are cell batches (plus ghosted ones for the MPI case).
> So in your code, you would fill the rho and mu fields by a cell loop (for
> MPI, please run it for "n_cells = this->data->n_macro_cells() +
> this->data->n_ghost_cell_batches()". Then you can access that data both on
> cells and faces, i.e., with FEEvaluation and FEFaceEvaluation.
>
> The implementation resolves the `interior_cells` and `exterior_cells`
> fields for you.
>
> Retrieving the `user_index` is a bit more involved as we do not provide
> access to a face iterator.
>
> You could do it approximately as follows (sitting on a face like in your
> code below):
>
> const unsigned int interior_cell = face2cell_info.cells_interior[element];
> const Triangulation<dim>::cell_iterator my_cell =
> this->data->get_cell_iterator(interior_cell /
> VectorizedArray<Number>::n_array_elements,
> interior_cell %
> VectorizedArray<Number>::n_array_elements);
> const Triangulation<dim>::face_iterator face_it =
> my_cell->face(face2cell_info.interior_face_no);
>
> Can you try if that works? If yes, we should add some access functions to
> MatrixFree. Let us know if you are interested in helping with a pull
> request.
>
Ok, it works (at least on active level). I tested that by comparing
location of quadrature points (code below). But for implementing
get_face_iterator I think that it would be safer to do it using
cell_level_index (eg. numbering of cell may change )
FEFaceEvaluation<dim, degree_u, degree_u+1+0 , dim, Number>
phi_inner(*(this->data),
true);
FEFaceValues<dim> fe_face_values(this->data->get_dof_handler().get_fe(),
this->data->get_face_quadrature(),
UpdateFlags::update_quadrature_points );
const unsigned int n_faces=this->data->n_inner_face_batches();
mu_inner_faces.resize(n_faces);
for (unsigned int face=0; face<n_faces; ++face){
phi_inner.reinit(face);
const internal::MatrixFreeFunctions::FaceToCellTopology
< VectorizedArray< Number >::n_array_elements >
face2cell_info=this->data->get_face_info(face);
for(unsigned int element=0;
element<this->data->n_active_entries_per_face_batch(face) ;
++element){
const unsigned int interior_cell = face2cell_info.cells_interior[element];
// const unsigned int exterior_cell =
face2cell_info.cells_exterior[element];
const typename Triangulation<dim>::cell_iterator my_cell =
this->data->get_cell_iterator(interior_cell /
VectorizedArray<Number>::n_array_elements,
interior_cell %
VectorizedArray<Number>::n_array_elements);
const typename Triangulation<dim>::face_iterator face_it =
my_cell->face(face2cell_info.interior_face_no);
fe_face_values.reinit(my_cell, face2cell_info.interior_face_no);
std::vector<Point< dim> > iterator_points =
fe_face_values.get_quadrature_points ();
for(unsigned int point_number=1; point_number<iterator_points.size();
++point_number){
Point< dim, VectorizedArray< Number > > point_face
=phi_inner.quadrature_point (point_number);
double difference =0;
for(unsigned int d=0; d<dim; ++d){
difference+=
std::fabs(iterator_points[point_number](d)-point_face(d)[element]);
}
if (difference>1e-8)
std::cout<<"Diffrence: " <<difference<< " in quadrature points on
faces!"<<std::endl;
}
mu_inner_faces[face][element]=1;
}
}
--
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/13c24a74-caad-4c04-9b9c-d68e7b874d4c%40googlegroups.com.
