Skip to content

Fe br #8

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
bangerth wants to merge 10 commits into
base: master
Choose a base branch
from
Open

Fe br #8

Show file tree
Hide file tree
Changes from 8 commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
c3ca881
Add functionality to support multiple mappings per element
GrahamBenHarper Apr 24, 2019
4873bf9
Generalize mapping types per DoF
GrahamBenHarper Apr 25, 2019
d83be4b
Correct rt/nedelec sign change and update documentation and indentation
GrahamBenHarper Apr 26, 2019
4dcecfd
Update FE_BernardiRaugel mapping
GrahamBenHarper Apr 26, 2019
1f53354
Change fe_bernardi_raugel bubble function maps to mapping_piola
GrahamBenHarper Apr 27, 2019
19cb379
Allocate for multiple mappings in fe_poly_tensor::get_data
GrahamBenHarper May 24, 2019
d180b1a
Adjust fe_poly_tensor::get_data to allocate memory more efficiently
GrahamBenHarper May 26, 2019
760b276
Add override to Bernardi-Raugel fill_fe_values but still call base fu…
GrahamBenHarper Jun 18, 2019
1a9ba9f
Polish FE_PolyTensor::get_data and other commands so they are more clear
GrahamBenHarper Jun 26, 2019
c14d37f
Run indent script
GrahamBenHarper Jun 27, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 20 additions & 0 deletions include/deal.II/fe/fe_bernardi_raugel.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,6 +102,26 @@ class FE_BernardiRaugel
*/
void
initialize_support_points();

/**
* Fill the values Bernardi-Raugel polynomials, but replace the normal
* vector from the PolynomialsBernardiRaugel class with the physical
* normal vector from @p cell.
*/
virtual void
fill_fe_values(
const typename Triangulation<dim, dim>::cell_iterator &cell,
const CellSimilarity::Similarity cell_similarity,
const Quadrature<dim> & quadrature,
const Mapping<dim, dim> & mapping,
const typename Mapping<dim, dim>::InternalDataBase &mapping_internal,
const dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
dim>
& mapping_data,
const typename FiniteElement<dim, dim>::InternalDataBase &fe_internal,
dealii::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
dim>
&output_data) const override;
};

DEAL_II_NAMESPACE_CLOSE
Expand Down
2 changes: 1 addition & 1 deletion include/deal.II/fe/fe_dg_vector.templates.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@ FE_DGVector<PolynomialType, dim, spacedim>::FE_DGVector(const unsigned int deg,
std::vector<ComponentMask>(PolynomialType::compute_n_pols(deg),
ComponentMask(dim, true)))
{
this->mapping_type = map;
this->mapping_type = {map};
const unsigned int polynomial_degree = this->tensor_degree();

QGauss<dim> quadrature(polynomial_degree + 1);
Expand Down
38 changes: 30 additions & 8 deletions include/deal.II/fe/fe_poly_tensor.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -207,9 +207,12 @@ class FE_PolyTensor : public FiniteElement<dim, spacedim>
protected:
/**
* The mapping type to be used to map shape functions from the reference
* cell to the mesh cell.
* cell to the mesh cell. If this vector is length one, the same mapping
* will be applied to all shape functions. If the vector size is equal to
* the finite element dofs per cell, then each shape function will be mapped
* according to the corresponding entry in the vector.
*/
MappingType mapping_type;
std::vector<MappingType> mapping_type;


/* NOTE: The following function has its definition inlined into the class
Expand Down Expand Up @@ -245,11 +248,33 @@ class FE_PolyTensor : public FiniteElement<dim, spacedim>
// initialize fields only if really
// necessary. otherwise, don't
// allocate memory

bool update_transformed_shape_values = false;
bool update_transformed_shape_grads = false;
bool update_transformed_shape_hessian_tensors = false;

for (unsigned int i = 0; i < this->mapping_type.size(); ++i)
{
MappingType mapping_type = this->mapping_type[i];

if (mapping_type != mapping_none)
update_transformed_shape_values = true;
Copy link
Owner Author

@bangerth bangerth Jun 21, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In essence, what you're doing here is determine whether any of the elements of the vector are different from mapping_none. (And similarly in the next few lines.) You can write this more concisely in the following way (see https://en.cppreference.com/w/cpp/algorithm/find):

  const book update_transformed_shape_values = (std::find_if_not (this->mapping_type.begin(),
                                                                                                 this->mapping_type.begin(),
                                                                                                 [](const MappingType t) { return t != mapping_none; })
                                                                                   != this->mapping_type.end());

Copy link

@GrahamBenHarper GrahamBenHarper Jun 26, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This command is a mouthful. Let me see... so this command locates any entry between begin or end satisfies the inline function (what's the C++ word for that?) that returns true if the mapping type is not mapping_none. I assume the default return is .end() if nothing is found?

 const bool update_transformed_shape_values = (std::find_if (this->mapping_type.begin(),
                                                             this->mapping_type.end(),
                                                             [](const MappingType t) { return t != mapping_none; })
                                                                                   != this->mapping_type.end());

Should that be .end() on the second line? Also, if we use return t != mapping_none, shouldn't we use find_if instead of find_if_not?

And if I were to do this for transformed_shape_grads, would it look like this?

 const bool update_transformed_shape_values = (std::find_if (this->mapping_type.begin(),
                                                             this->mapping_type.end(),
                                                             [](const MappingType t) { return (t == mapping_raviart_thomas ||
                                                                                               t == mapping_piola ||
                                                                                               t == mapping_nedelec ||
                                                                                               t == mapping_contravariant); })
                                                                                   != this->mapping_type.end());

If so, then I think I understand how this works and I can implement it.

Copy link

@GrahamBenHarper GrahamBenHarper Jun 26, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh is that called a "lambda expression"? This is something I haven't explored much of before.

Copy link
Owner Author

@bangerth bangerth Jul 1, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes:

  • That's a lambda function -- a function without a name defined in the place where you need it.
  • And yes, correct, it should have been .end() in the second line.
  • And yes, if find_if can't find an entry that satisfies the lambda function, then it returns .end().

Copy link
Owner Author

@bangerth bangerth Jul 1, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

https://en.wikipedia.org/wiki/Anonymous_function#C%2B%2B_%28since_C%2B%2B11%29


if ((mapping_type == mapping_raviart_thomas) ||
(mapping_type == mapping_piola) ||
(mapping_type == mapping_nedelec) ||
(mapping_type == mapping_contravariant))
update_transformed_shape_grads = true;
Copy link
Owner Author

@bangerth bangerth Jun 21, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This would work with find_if (...).


if (mapping_type != mapping_none)
update_transformed_shape_hessian_tensors = true;
Copy link
Owner Author

@bangerth bangerth Jun 21, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is the exact same condition as you had above for update_transformed_shape_values, so you will get the exact same answer. Is this on purpose?

Copy link

@GrahamBenHarper GrahamBenHarper Jun 26, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah it was not on purpose. Seeing the original command looked like

            if (mapping_type != mapping_none)
              data->untransformed_shape_hessian_tensors.resize(n_q_points);

I can just absorb that boolean as being equal to the update_transformed_shape_values boolean since both only check if the mapping is not mapping_none.

}

if (update_flags & update_values)
{
values.resize(this->dofs_per_cell);
data->shape_values.reinit(this->dofs_per_cell, n_q_points);
if (mapping_type != mapping_none)
if (update_transformed_shape_values)
data->transformed_shape_values.resize(n_q_points);
}

Expand All @@ -259,10 +284,7 @@ class FE_PolyTensor : public FiniteElement<dim, spacedim>
data->shape_grads.reinit(this->dofs_per_cell, n_q_points);
data->transformed_shape_grads.resize(n_q_points);

if ((mapping_type == mapping_raviart_thomas) ||
(mapping_type == mapping_piola) ||
(mapping_type == mapping_nedelec) ||
(mapping_type == mapping_contravariant))
if (update_transformed_shape_grads)
data->untransformed_shape_grads.resize(n_q_points);
}

Expand All @@ -271,7 +293,7 @@ class FE_PolyTensor : public FiniteElement<dim, spacedim>
grad_grads.resize(this->dofs_per_cell);
data->shape_grad_grads.reinit(this->dofs_per_cell, n_q_points);
data->transformed_shape_hessians.resize(n_q_points);
if (mapping_type != mapping_none)
if (update_transformed_shape_hessian_tensors)
data->untransformed_shape_hessian_tensors.resize(n_q_points);
}

Expand Down
2 changes: 1 addition & 1 deletion source/fe/fe_abf.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,7 @@ FE_ABF<dim>::FE_ABF(const unsigned int deg)
Assert(dim >= 2, ExcImpossibleInDim(dim));
const unsigned int n_dofs = this->dofs_per_cell;

this->mapping_type = mapping_raviart_thomas;
this->mapping_type = {mapping_raviart_thomas};
// First, initialize the
// generalized support points and
// quadrature weights, since they
Expand Down
2 changes: 1 addition & 1 deletion source/fe/fe_bdm.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,7 @@ FE_BDM<dim>::FE_BDM(const unsigned int deg)

const unsigned int n_dofs = this->dofs_per_cell;

this->mapping_type = mapping_bdm;
this->mapping_type = {mapping_bdm};
// These must be done first, since
// they change the evaluation of
// basis functions
Expand Down
48 changes: 47 additions & 1 deletion source/fe/fe_bernardi_raugel.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,25 @@ FE_BernardiRaugel<dim>::FE_BernardiRaugel(const unsigned int p)

// const unsigned int n_dofs = this->dofs_per_cell;

this->mapping_type = mapping_none;
if(dim==2)
this->mapping_type = {mapping_none, mapping_none,
mapping_none, mapping_none,
mapping_none, mapping_none,
mapping_none, mapping_none,
mapping_piola, mapping_piola,
mapping_piola, mapping_piola};
Copy link
Owner Author

@bangerth bangerth Jun 21, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I believe you'll have to update that, right?

Copy link

@GrahamBenHarper GrahamBenHarper Jun 26, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I will revert it to {mapping_none} once the fill_fe_values overload is finished.

else if(dim==3)
this->mapping_type = {mapping_none, mapping_none, mapping_none,
mapping_none, mapping_none, mapping_none,
mapping_none, mapping_none, mapping_none,
mapping_none, mapping_none, mapping_none,
mapping_none, mapping_none, mapping_none,
mapping_none, mapping_none, mapping_none,
mapping_none, mapping_none, mapping_none,
mapping_none, mapping_none, mapping_none,
mapping_piola, mapping_piola,
mapping_piola, mapping_piola,
mapping_piola, mapping_piola};
Copy link
Owner Author

@bangerth bangerth Jun 21, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

add the following:

  else
    Assert (false, ExcNotImplemented());

to make sure that if the code is ever used for anything other than 2d and 3d, you will be notified that that case has not been implemented.

Copy link

@GrahamBenHarper GrahamBenHarper Jun 26, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

On line 53 I have Assert(dim == 2 || dim == 3, ExcImpossibleInDim(dim));

Should I still add this one as well?

Copy link
Owner Author

@bangerth bangerth Jul 1, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I think that would be useful. The assertion at the top is an expression of the general idea that this function isn't implemented for certain cases, whereas adding the assertion here would point to a specific place where could would need to be augmented if one wanted to extend the function.

Copy link

@GrahamBenHarper GrahamBenHarper Jul 4, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay got it!

// These must be done first, since
// they change the evaluation of
// basis functions
Expand Down Expand Up @@ -178,6 +196,34 @@ FE_BernardiRaugel<dim>::initialize_support_points()
}
}



template <int dim>
void
FE_BernardiRaugel<dim>::fill_fe_values(
const typename Triangulation<dim, dim>::cell_iterator &cell,
const CellSimilarity::Similarity cell_similarity,
const Quadrature<dim> & quadrature,
const Mapping<dim, dim> & mapping,
const typename Mapping<dim, dim>::InternalDataBase &mapping_internal,
const dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
dim>
& mapping_data,
const typename FiniteElement<dim, dim>::InternalDataBase &fe_internal,
dealii::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
dim>
&output_data) const
{
FE_PolyTensor<PolynomialsBernardiRaugel<dim>, dim, dim>::fill_fe_values(cell,
cell_similarity,
quadrature,
mapping,
mapping_internal,
mapping_data,
fe_internal,
output_data);
}
Copy link
Owner Author

@bangerth bangerth Jun 21, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, that's what I had in mind. Now you'll just have to overwrite what this function has computed for the four last shape functions.

Copy link

@GrahamBenHarper GrahamBenHarper Jun 26, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes I like that idea. I stopped here because I need to read what kind of data or structure is sitting in output_data so I know what I need to do to modify it properly.


template class FE_BernardiRaugel<1>;
template class FE_BernardiRaugel<2>;
template class FE_BernardiRaugel<3>;
Expand Down
6 changes: 3 additions & 3 deletions source/fe/fe_dg_vector.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@ DEAL_II_NAMESPACE_OPEN

template <int dim, int spacedim>
FE_DGNedelec<dim, spacedim>::FE_DGNedelec(const unsigned int p)
: FE_DGVector<PolynomialsNedelec<dim>, dim, spacedim>(p, mapping_nedelec)
: FE_DGVector<PolynomialsNedelec<dim>, dim, spacedim>(p, {mapping_nedelec})
{}


Expand All @@ -52,7 +52,7 @@ template <int dim, int spacedim>
FE_DGRaviartThomas<dim, spacedim>::FE_DGRaviartThomas(const unsigned int p)
: FE_DGVector<PolynomialsRaviartThomas<dim>, dim, spacedim>(
p,
mapping_raviart_thomas)
{mapping_raviart_thomas})
{}


Expand All @@ -77,7 +77,7 @@ FE_DGRaviartThomas<dim, spacedim>::get_name() const

template <int dim, int spacedim>
FE_DGBDM<dim, spacedim>::FE_DGBDM(const unsigned int p)
: FE_DGVector<PolynomialsBDM<dim>, dim, spacedim>(p, mapping_bdm)
: FE_DGVector<PolynomialsBDM<dim>, dim, spacedim>(p, {mapping_bdm})
{}


Expand Down
2 changes: 1 addition & 1 deletion source/fe/fe_nedelec.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,7 @@ FE_Nedelec<dim>::FE_Nedelec(const unsigned int order)

const unsigned int n_dofs = this->dofs_per_cell;

this->mapping_type = mapping_nedelec;
this->mapping_type = {mapping_nedelec};
// First, initialize the
// generalized support points and
// quadrature weights, since they
Expand Down
Loading