[cig-commits] commit 1866 by bangerth to /var/svn/dealii/aspect

[dealii.demon at gmail.com]

Revision 1866

Patch suggested by an email by Thomas Geenen: Use a more reasonable number of quadrature points when assembling the linear systems for temperature or compositional fields.

U   trunk/aspect/doc/modules/changes.h
U   trunk/aspect/source/simulator/assembly.cc


http://www.dealii.org/websvn/revision.php?repname=Aspect+Repository&path=%2F&rev=1866&peg=1866

Diff:
Modified: trunk/aspect/doc/modules/changes.h
===================================================================
--- trunk/aspect/doc/modules/changes.h	2013-08-30 16:47:55 UTC (rev 1865)
+++ trunk/aspect/doc/modules/changes.h	2013-08-30 22:33:50 UTC (rev 1866)
@@ -8,6 +8,12 @@
 </p>
 
 <ol>
+  <li>Fixed: The assembly of the temperature and compositional linear systems
+  in each time step used an unnecessarily large number of quadrature points.
+  This is now fixed.
+  <br>
+  (Thomas Geenen, 2013/08/30)
+
   <li>Extended: The ability to compute with tracers has been completely
   overhauled. In particular, there is now a cookbook in the manual that
   describes how to use them.

Modified: trunk/aspect/source/simulator/assembly.cc
===================================================================
--- trunk/aspect/source/simulator/assembly.cc	2013-08-30 16:47:55 UTC (rev 1865)
+++ trunk/aspect/source/simulator/assembly.cc	2013-08-30 22:33:50 UTC (rev 1866)
@@ -1382,8 +1382,24 @@
                           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.is_temperature()
+                                            ?
+                                            parameters.temperature_degree
+                                            :
+                                            parameters.composition_degree)
+                                           +
+                                           (parameters.stokes_velocity_degree+1)/2),
                                parameters.n_compositional_fields),
          internal::Assembly::CopyData::
          AdvectionSystem<dim> (finite_element));