FE_PartitionOfUnity

[Denis Davydov]

Dear all, 

I have been developing a partition of unity FE for deal.ii, i.e.

U(x) = \sum_i N_i(x) U_i + \sum_j N_j(x) \sum_k F_k(x) U_{jk}

where N_i and N_j are possibly different FE bases (including a mapping from the reference cell) 

and F_k(x) are the enrichment functions in real space (1/r, 1/r^2, Heaviside, etc); 

U_i and U_{jk} are DoFs.

I pretty much follow the yet unpublished step-47 https://www.dealii.org/developer/doxygen/deal.II/step_47.html by exploiting FESystem.

So far for simplicity I derive from the FESystem and make everything related to FEValues and FEFaceValues work fine

by implementing fill_XYZ_values which call FESystem's ones and do some extra work.

I should be able to assemble weak forms now.

My current questions are:

1. how does a FE know that it's a vector-valued? I surely want it to be a scalar valued, but

   having derived from FESystem makes DataOut split output into the base DoFs and those which are enriched using a  given function in space.

   Consequently I also can't run get_function_values(...,std::vector<double>) due to the failing assert inside.

2. I tried outputting shape functions and clearly see the multiplication by  enrichment function is missing. 

   Which functions are actually called inside DataOut? By trial-and-error I can say that it is not

   FiniteElement::shape_value() and neither FiniteElement::fill_fe_values().

   

3. I clearly can't implement FiniteElement::shape_value(const unsigned int i, const Point<dim>& p) for arbitrary i 

   as p is a point on the reference element, whereas for partition of unity one evaluates a given

   function in real space. Could that be a dealbreaker? Is this function used anywhere?

4. Can I still have n_base_elements() >= 2  , n_components() = 1 , n_blocks() >= 2? That's effectively what it is - 

   it's build from (at least) 2 different FEs, where the second one is multiplied by a given function in space.

   Or should I just give up any re-ordering possibility according to the underlying FEs (no idea if this could be useful)?

So far I am leaning towards aggregation (keep FESystem inside in the private block and use what's relevant within the class), 

as deriving from FESystem leads to a vector-valued FE, which is not what I want. 

Any advices are much appreciated,

Regards,

Denis.

[Wolfgang Bangerth]

> 1. how does a FE know that it's a vector-valued? I surely want it to be
> a scalar valued, but
> having derived from FESystem

Yes, this won't work. FESystem counts the base elements and
multiplicities and tells the base class FiniteElement how many vector
components that amounts to. I think you'll have to copy the logic of
FESystem to handle base elements, but change what you communicate to the
base class.

Or you could make a FESystem a member variable of your new class. You'd
then call the FiniteElement base constructor yourself, but in all other
functions essentially refer to the member variable.

> 2. I tried outputting shape functions and clearly see the multiplication
> by enrichment function is missing.
> Which functions are actually called inside DataOut? By
> trial-and-error I can say that it is not
> FiniteElement::shape_value() and neither
> FiniteElement::fill_fe_values().

I think DataOut goes through FEValues which itself calls
FiniteElement::fill_fe_values() on the relevant element in use.

> 3. I clearly can't implement FiniteElement::shape_value(const unsigned
> int i, const Point<dim>& p) for arbitrary i
> as p is a point on the reference element, whereas for partition of
> unity one evaluates a given
> function in real space. Could that be a dealbreaker? Is this
> function used anywhere?

That's one of the issues I'm currently looking through in #1198. I think
that for the most part, these functions that compute values of shape
functions are not widely used. Certainly not for integration where
everything goes through FEValues which itself goes through
fill_fe_values. But there are functions in FETools that compute, for
example, interpolation matrices. These may simply not work if you don't
implement the functions in question.

> 4. Can I still have n_base_elements() >= 2 , n_components() = 1 ,
> n_blocks() >= 2? That's effectively what it is -

Good question. That's worth trying.

> So far I am leaning towards aggregation (keep FESystem inside in the
> private block and use what's relevant within the class),
> as deriving from FESystem leads to a vector-valued FE, which is not what
> I want.

Yes, this sounds reasonable.

Best
W>

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@math.tamu.edu
www: http://www.math.tamu.edu/~bangerth/

[Denis Davydov]

Thanks Wolfgang for the prompt reply.

On Wednesday, August 12, 2015 at 8:16:35 PM UTC+2, bangerth wrote:

Or you could make a FESystem a member variable of your new class. You'd
then call the FiniteElement base constructor yourself, but in all other
functions essentially refer to the member variable.

could you please elaborate what each argument in constructor does?

The documentation states "Constructor", which is not very descriptive ;-)

I am fine with the FiniteElementData class, but two others are not obvious,

namely :

const std::vector< bool > & restriction_is_additive_flags,

const std::vector< ComponentMask > & nonzero_components . 

What should their size be?

Or if there are places where the two other arguments are discussed in details, 

perhaps one could add a link in documentation of the constructor.

 

> So far I am leaning towards aggregation (keep FESystem inside in the
> private block and use what's relevant within the class),
> as deriving from FESystem leads to a vector-valued FE, which is not what
> I want.

Yes, this sounds reasonable.

Great. will go this route.

Kind regards,

Denis.

[Wolfgang Bangerth]

> Or you could make a FESystem a member variable of your new class. You'd
> then call the FiniteElement base constructor yourself, but in all other
> functions essentially refer to the member variable.
>
>
> could you please elaborate what each argument in constructor does?
> The documentation states "Constructor", which is not very descriptive ;-)

Yes, these old parts of the library are not well documented :-(

> I am fine with the FiniteElementData class, but two others are not obvious,
> namely :
> const std::vector< bool > & restriction_is_additive_flags,

The best I can currently provide is probably in this link to
get_restriction_matrix():

https://www.dealii.org/developer/doxygen/deal.II/classFiniteElement.html#a75f248154c84e88e2fd1a8946647d262

and this member variable documentation (where the constructor stores the
argument):

https://www.dealii.org/developer/doxygen/deal.II/classFiniteElement.html#aa8f4833a318127b42d2dc806bffe1c2c

The array needs to have size dofs_per_cell, and you probably simply need to
take for each of your base functions the corresponding value of your base
element (in the same way as FESystem must likely do).

> const std::vector< ComponentMask > & nonzero_components .

This denotes for each shape function in which vector components it is nonzero.
In your case, you have as many vector components as your base elements have.
Again, you should just take the component masks from the shape functions of
your base elements. (FESystem will here take the component masks from the base
elements, but then expand them to the size of the whole FESystem's number of
vector components.)

Best
W.

[Denis Davydov]

Thanks, Wolfgang. 

I think i made the switch to aggregation without breaking fill_fe_(face/subface)_values i had working before.

Will go on with h- and hp-related functions.

Regards,

Denis.

[Wolfgang Bangerth]

> could you please elaborate what each argument in constructor does?
> The documentation states "Constructor", which is not very descriptive ;-)
> I am fine with the FiniteElementData class, but two others are not obvious,
> namely :
> const std::vector< bool > & restriction_is_additive_flags,
> const std::vector< ComponentMask > & nonzero_components .
> What should their size be?
> Or if there are places where the two other arguments are discussed in details,
> perhaps one could add a link in documentation of the constructor.

There was nothing, but you will be interested in this pull request now:
https://github.com/dealii/dealii/pull/1398

Best
W.

[Denis Davydov]

I think i have issues with CellSimilarity::Similarity and update_flags. 

If for the case of fill_fe_values i use flags from 

  const UpdateFlags flags(fe_data.current_update_flags());

then on the second cell update_values is not in current_update_flags() and

so i can't do something like: if (flags & update_values) 

If I use instead

  const UpdateFlags flags(fe_data.update_flags);

then, supposedly due to CellSimilarity, the shape values provided in 

internal::FEValues::FiniteElementRelatedData<dim,spacedim> &output_data

on the 2nd cell are those cached from the previous cell (not re-calculated). 

In other words, say for Trapezoidal rule and linear FEs on the first cell i have

  element_dof=7 N=1 f=0.326922 N*f=0.326922

then on the next cell the shape value at the same element-quadrature point q is the result of previous calculation:

  element_dof =7 N=0.326922 f=0.243117 N*f=0.0794802

whereas what i want to have is the "vanilla" shape function cached, i.e. N=1 above. bummer...

I suppose a workaround is to reset CellSimilarity inside fill_fe_values (and others) before I call

fe_system.fill_fe_values(mapping,

                           cell,

                           quadrature,

                           mapping_internal,

                           fe_internal,

                           mapping_data,

                           output_data,

                           cell_similarity);

and perhaps adding update_values and upate_gradients to flags inside update_each() 

to make sure fe_system indeed re-calculate shape values and gradients for each cell

so that I can safely post-multiply those using the enrichment function.

Any better ways to do it?

Regards,

Denis.

 

[Denis Davydov]

On Thursday, August 20, 2015 at 6:24:43 PM UTC+2, bangerth wrote:

There was nothing, but you will be interested in this pull request now:
   https://github.com/dealii/dealii/pull/1398

Thanks, Wolfgang! This makes things much more clear w.r.t. FEValues class.

Regards,

Denis. 

[Denis Davydov]

On Thursday, August 20, 2015 at 6:26:01 PM UTC+2, Denis Davydov wrote:

Any better ways to do it?

This didn't work, hm...

[Denis Davydov]

On Thursday, August 20, 2015 at 8:07:00 PM UTC+2, Denis Davydov wrote:

This didn't work, hm...

Nevermind, Solved it. 

Had to pull data from base_fe_data and reset whatever was cached from the previous cell,

effectively copy-paste with minor modifications of the bottom part of FESystem::compute_fill_one_base().

h-refinement seems to be working. 

Will look at hp now...

Regards,

Denis.

[Denis Davydov]

A question of hp:

for face domination something like

Q(2) vs Q(4)  leads to Q(2) dominating Q(4). All good.

something like 

Q(2) vs Q(4) + enriched Q(1)  is still good, 

i constrain some DoFs (base and enriched) of the second element to make the whole thing continuous.

However, 

Q(4) vs Q(2) + enriched Q(1) is now a problem. 

On one hand, i want to constrain dofs of Q4 to match the Q2. 

On another hand, i also want to zero some of the enriched dofs to make the field continuos.

I wonder if I can handle this situation at all somehow? 

This situation could easily happen when enriched elements are fine and do not need higher polynomial degree, 

whereas some neighbouring non-enriched elements may benefit from pumping their degree.

The only alternative I see so far is to implement constrains on enriched DoFs in a separate function, 

which loops over all enriched elements and if there is non-enriched element sitting nearby  

constrain some DoFs to zero. That's a bit ugly but should work in principle. I still hope to avoid it...

Regards,

Denis.

[Wolfgang Bangerth]

Denis,

> for face domination something like
> Q(2) vs Q(4) leads to Q(2) dominating Q(4). All good.
>
> something like
> Q(2) vs Q(4) + enriched Q(1) is still good,
> i constrain some DoFs (base and enriched) of the second element to make the
> whole thing continuous.
>
> However,
> Q(4) vs Q(2) + enriched Q(1) is now a problem.
> On one hand, i want to constrain dofs of Q4 to match the Q2.
> On another hand, i also want to zero some of the enriched dofs to make the
> field continuos.

Correct. I don't quite recall how we do this in the library (if at all), but
it's essentially the same case as when you have a large Q4 element meeting two
small Q2 elements along a refined edge -- neither the restriction of the Q4
nor the restriction of the two Q2s dominate each other. You need to find a
third space that dominates both, which is the Q2 on the large edge. There is a
section about this in the hp paper by myself and Kayser-Herold on exactly this
case. We discuss hanging nodes, but the same arguments will be true for the
case you consider.

> I wonder if I can handle this situation at all somehow?

As long as the pure Q2 is in your hp::FECollection, this could be resolved.
I'm not sure the case is implemented (see above), but I don't see an
overwhelming obstacle to making this happen.

The beginning of a related discussion (for the opposite direction: not
included but including spaces) is here, by the way:
https://github.com/dealii/dealii/issues/313

Best
Wolfgang

[Denis Davydov]

On Saturday, August 22, 2015 at 5:27:00 AM UTC+2, bangerth wrote:

As long as the pure Q2 is in your hp::FECollection, this could be resolved.
I'm not sure the case is implemented (see above), but I don't see an
overwhelming obstacle to making this happen.

for a record, here is a PR with a couple of steps in that direction: https://github.com/dealii/dealii/pull/1424

Regards,

Denis. 

[Denis Davydov]

After those PRs are out of the way, 

I think there is another challenging thing:  rescaling of enriched shape functions to have better condition number.

Say I have f(x)=exp(-x) as enrichment in 1D on top of linear shape functions N(x). 

For elements close to 0 the enriched shape functions ( N(x) f(x) ) would be not too bad, 

but, say, at a distance 100 the local mass matrix will be terrible as N(x) exp(-x)

will be close to zero on that element...

For a moment let's assume i just scale N(x) exp(-x) so that it satisfies delta property.

That's certainly doable, in places where I have something like:

   output_data.shape_values[enriched_to_complete[d]][q] *= enrichment->value(mapping_data.quadrature_points[q]);

I can keep a map  {global_dof} -> scalar multiplier and the multiply by a scalar.

I would need to pre-calculate it for a given DoFHandler before doing anything, though...

Interesting part are h-constraints...
Say for h-case only and (un-enriched) linear elements I have

c2 = (c0+c1)/2 = [0.5, 0.5] [c0, c1].

the same constraint hold if I do enrichment BUT without any normalisation.

Say I decide to normalise each enriched DoF and introduce

d2 = a2 * c2

d0 = a0 * c0

d1 = a1 * c1

with non-zero a_i.

I suppose I can loop over constraint matrix and correct its entries for enriched DoFs such that:

d2 = [0.5 * a2/a0,  0.5 * a2/a1] [d0, d1]

which should make things continuous even with local renormalisation.

That's a fife-minute quick though, but so far I think it should work... at least for h-refinement with linear elements.

Also given the fact that those enriched DoFs live on a separate "level" (1st component in the underlying FESystem),

then they are constrained only between each other. So in constraint matrix I need to touch only those entries which

constrain enriched DoFs and safely skip everything else.

What do you think?

Regards,

Denis.

[Denis Davydov]

Seems to be working (see attached JPGs for shape functions on [0,1]^2 corresponding to bilinear shape functions with exp(-|x|) enrichment).

What I currently do is quite ugly - I copy a ConstraintMatrix and recreate it based on a stored 

std::map of  global POU dofs -> normalisation value.

Ideally, that should be inside DoFTools::make_hanging_node_constraints(), but it will not fit the current logic

as functions like get_face_interpolation_matrix() are agnostic to the face (and thus its location in real space) 

on which they are to be used.

A workaround is to create some function in the constraint matrix which will re-scale constraints as needed.

I would happily do it in a separate PR, if you think that this is a good way to do it.

That would be something like:

ConstraintMatrix.rescale(const std::map<types::global_dof_index, double> scaling_factors);

One can also see it as a multiplication with a diagonal matrix, which has quite few elements not equal to 1.

(In real-life most DoFs will be untouched that is not scaled).

I would then be using this function after:

1) I calculate scaling factors for each DoF ( static function of the my FiniteElement class which receives DoFHandler, workout support points for POU DoFs and get needed scaling factor for given enrichment function)

2) DoFTools::make_hanging_node_constraints  (dof_handler, constraints) is called.

Regards,

Denis.

[Denis Davydov]

an update to the topic:

I decided to avoid re-normalisation. Otherwise, the class seems to be working fine. I successfully solve several test problems with enrichment. 

I plan to clean things up and submit a PR in February.

Regards,

Denis.

[Denis Davydov]

In case someone will come across this topic, 

FE_Enriched class has been added to the deal.II library, see https://www.dealii.org/developer/doxygen/deal.II/classFE__Enriched.html

Regards,

Denis.