[aspect-devel] matrix assembly

[Wolfgang Bangerth]

> The reasoning may have been that we have to assemble a term of the form
>    u . grad phi_i * phi_j
> with a total polynomial degree of
>    stokes_degree + 2 * temp_degree - 1
> which would suggest a Gauss formula of degree
>    temp_degree + stokes_degree/2
> rounded up if necessary. Does this sound reasonable? It would boil down to
>    temp_degree + 1
> in the most common case.

This would be the attached patch. Makes sense?
W.

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

-------------- next part --------------
Index: source/simulator/assembly.cc
===================================================================
--- source/simulator/assembly.cc	(revision 1847)
+++ source/simulator/assembly.cc	(working copy)
@@ -1382,8 +1382,25 @@
                           copy_local_to_global_advection_system,
                           this,
                           std_cxx1x::_1),
+
+                          // we have to assemble the term u.grad phi_i * phi_j, which is
+                          // of total polynomial degree
+                          //   stokes_deg + 2*temp_deg -1
+                          // (or similar for comp_deg). this suggests using a Gauss
+                          // quadrature formula of order
+                          //   temp_deg + stokes_deg/2
+                          // rounded up. do so. (note that x/2 rounded up
+                          // equals (x+1)/2 using integer division.)
          internal::Assembly::Scratch::
-         AdvectionSystem<dim> (finite_element, mapping, QGauss<dim>(parameters.composition_degree+2),
+         AdvectionSystem<dim> (finite_element, mapping,
+                               QGauss<dim>((temperature_or_composition ==
+                                            TemperatureOrComposition::temperature_field
+                                            ?
+                                            parameters.temperature_degree
+                                            :
+                                            parameters.composition_degree)
+                                           +
+                                           (parameters.stokes_velocity_degree+1)/2),
                                parameters.n_compositional_fields),
          internal::Assembly::CopyData::
          AdvectionSystem<dim> (finite_element));