Add a push_back member function to Quadrature<dim>.

[simonsticko]

Immersed quadrature rules are complicated to construct. It is often
done in several steps over different parts of a cell. The
implementation of the algorithm generating immersed quadratures will be
simplified if we can pass around a reference to a Quadrature and
push back quadrature points on it consecutively.

[simonsticko]

I have some code that generates immersed quadrature rules from a level set function. This is one step towards getting it into dealii.

[drwells]

I am not sure that this is the right approach: it implies that the Quadrature object is in some undetermined state before and after we call this function where it is not quite a quadrature rule. What is the advantage to using push_back instead of just accumulating into a pair of arrays and then calling Quadrature::initialize?

[simonsticko]

It's just intended to make it easier to construct quadrature rules.

In X/cut/immersed-FEM one often wants to create a quadrature rule over part of a cell. Creating these is typically quite involved. See for example

http://www.sciencedirect.com/science/article/pii/S0045782516308696
http://epubs.siam.org/doi/10.1137/140966290

Often generation is done by some form of subdivision of the cell, where points are created in several steps. Having a push_back function makes the implementation easier, since we can then pass around a quadrature and add points to it in several steps.

One could of course pass around something like

template<int dim>
struct IncompleteQuadrature{
    std::vector<Point<dim>> points;
    std::vetor<double> weights;
}

and then create a Quadrature from this. However, this is essentially the same container as the Quadrature class, so I think it makes more sense to use the container that is already in dealii.

One could argue that allowing for a push_back function allows the user to modify (destroy) a quadrature rule that is already correct. However, whether or not the Quadrature should be modified or not can be handled through constness.

[bangerth]

I can agree with @drwells 's objections. On the other hand, there is relatively little harm to adding the function. We almost always store const Quadrature objects for which the new function can not be called.

[simonsticko]

I didn't like my first commit message so I rewrote it.

[bangerth]

I think I'm ok with this patch, but I'm not going to overrule @drwells's objection. @drwells -- what do you think?

[drwells]

I thought about it some more and I still think there has to be a better way to solve the problem than putting this function in the base class. I would prefer to implement

template <int dim>
class ImmersedQuadrature : public Quadrature<dim>
{
public:
  void push_back(double weight, Point<dim> point);
};

so that one can assemble immersed quadratures point-by-point in a convenient way but not add an extra function to the base class.

I found another example of a problem created by adding this function:
ImmersedSurfaceQuadrature already has its own push_back that also requires three arguments instead of two (a point, a weight, and a normal vector). Hence, calling the base class' push_back will put objects in a broken state. We can avoid this issue by putting the two-argument version in a different class hierarchy.

[simonsticko]

I would be okey with the solution suggested by @drwells, but I would like to explain why I think this is an odd solution.

I see the Quadrature class just as a container for points an weights. It's not much different from std::vector<std::pair<Point<dim>,double>>, only more convenient to work with. From this point of view push_back makes sense because you clearly have to fill the container in some way.
Whether the Quadrature is sensible or not depends on what it contains. There are already constructors in Quadrature that put the weights to infinity. Surely that is an invalid state?

The suggested solution would be like having two related container classes, where you can only add elements to one of them:

template <class T>
class vector 
{
   ...
 // No push_back function here. 
};

template <class T>
ExtenableVector : public vector<T> {
public:
    void push_back(const T&);
};

But there is of course no need to do this because this can be handled through constness.

The mentioned issue with ImmersedSurfaceQuadrature exists regardless of this pull-request because you put ImmersedSurfaceQuadrature in a broken state through Quadrature<dim>::initialize.

I think the solution here is rather that ImmersedSurfaceQuadrature should not be derived from Quadrature. This also makes sense, because ImmersedSurfaceQuadrature<dim> is a way to integrate over a (dim-1)-dimensional set and not a dim-dimensional set. If I would pass an ImmersedSurfaceQuadrature to a function taking a Quadrature I would probably not get sensible results.

[drwells]

I thought about this some more and @simonsticko is right (it just took me awhile to understand why). We already have lots of ways to mutate Quadrature objects (and we even use them as collections of points with bogus weights in some contexts) so adding one more isn't a big deal.

[drwells]

/rebuild

[simonsticko]

Should I also set is_tensor_product_flag=false here?

[bangerth]

Yes. In fact it is this kind of issue that I think is so difficult to address if we make Quadrature mutable. My preference would have been if all members of Quadrature are at least conceptually const (they can't be in practice because we need operator=).

[simonsticko]

This problem already exist (if we really try to break something)
If we:

1. Create a tensor product quadrature.
2. Call Quadrature::initialize with something that isn't a tensor product.

we have a non-tensor product quadrature having is_tensor_product_flag=true. :(

The easy solution would be to just set is_tensor_product_flag=false in
Quadrature<dim>::initialize.

A more cumbersome solution (more along @drwells previous suggestion):

Since

1. push_back and initialize doesn't make much sense for a tensor product quadrature.
2. Not every quadrature "is a" tensor product.
   One could imagine the following

template<int dim>
class Quadrature : Subscriptor{
	
	// Same as before, but 
	// 1. initialize(..) removed.
	// 2. tensor_basis removed
	
	// 3. Changed to virtual
	virtual bool is_tensor_product(){ return false;} 
		
	// 4. Throws exception
	virtual const std::array<Quadrature<1>, dim> & get_tensor_basis() const;
};

template<int dim>
class TensorProductQuadrature : public Quadrature<dim> {
	// Constructors creating tensor products moved here.

        virtual bool is_tensor_product(){ return true;} 
	
	// Overloaded.
	const std::array<Quadrature<1>, dim> & get_tensor_basis() const override;
};

template<int dim>
class ExtendableQuadrature : Quadrature<dim>{
	
       // These only exist here
	void initialize(..);
	void push_back(..);
};

template<int dim>
class ImmersedSurfaceQuadrature : public Quadrature<dim>{

// This class no longer has push_back and initialize taking only points and
// weights. So we can't end up in a state with normals having a different size.
// Also, shorter implementation here since we can still derive from Quadrature.

};

[bangerth]

I like this class hierarchy. In fact, you don't even need the is_tensor_product() member function -- you can discover whether something is a tensor product quadrature using a dynamic_cast: i.e., the property is already encoded in the class hierarchy rather than via virtual member functions.

[simonsticko]

Since initialize is removed from Quadrature here, this is not really backwards compatible.
On the other hand, the number of users who create custom quadratures is probably small.

[bangerth]

Yes, I agree. Is that function used anywhere in the library?

[simonsticko]

I don't think so. If I remove it the library still compiles. I didn't run the whole test-suit but searching for
[q,Q][\w]+\.initialize\(
finds no occurrence of it.

[bangerth]

OK, that's a good indication that it's indeed unused. Thanks for trying!

[bangerth]

So what do we do with this patch?

[simonsticko]

If everyone agrees that the above solution is better I guess we just close this pull-request
and I'll open a new one with the solution above.