diff --git a/doc/modules/changes/20170513_egp-cig-reu b/doc/modules/changes/20170513_egp-cig-reu
new file mode 100644
index 00000000000..3504c952bdc
--- /dev/null
+++ b/doc/modules/changes/20170513_egp-cig-reu
@@ -0,0 +1,7 @@
+New particle interpolation scheme bilinear least squares
+using singular value decomposition algorithm. Currently
+only 2D models are supported. We chose a simple overshoot
+and undershoot correction scheme, mainly to truncate based
+on local cell maximum and minimum property value.
+
+(Elbridge Gerry Puckett, Ying He, Harsha Lokavarapu 2017/05/13)
\ No newline at end of file
diff --git a/include/aspect/particle/interpolator/bilinear_least_squares.h b/include/aspect/particle/interpolator/bilinear_least_squares.h
new file mode 100644
index 00000000000..a8b15083abf
--- /dev/null
+++ b/include/aspect/particle/interpolator/bilinear_least_squares.h
@@ -0,0 +1,62 @@
+/*
+ Copyright (C) 2017 by the authors of the ASPECT code.
+
+ This file is part of ASPECT.
+
+ ASPECT is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2, or (at your option)
+ any later version.
+
+ ASPECT is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with ASPECT; see the file doc/COPYING. If not see
+ .
+ */
+
+#ifndef _aspect_particle_interpolator_bilinear_least_squares_h
+#define _aspect_particle_interpolator_bilinear_least_squares_h
+
+#include
+#include
+
+namespace aspect
+{
+ namespace Particle
+ {
+ namespace Interpolator
+ {
+ /**
+ * Return the interpolated properties of all particles of the given cell using bilinear least squares method.
+ * Currently, only the two dimensional model is supported.
+ *
+ * @ingroup ParticleInterpolators
+ */
+ template
+ class BilinearLeastSquares : public Interface, public aspect::SimulatorAccess
+ {
+ public:
+ /**
+ * Return the cell-wise evaluated properties of the bilinear least squares function at the positions.
+ *
+ * @copydoc aspect::Particle::Interpolator::Interface::properties_at_points()
+ */
+ virtual
+ std::vector >
+ properties_at_points(const std::multimap > &particles,
+ const std::vector > &positions,
+ const ComponentMask &selected_properties,
+ const typename parallel::distributed::Triangulation::active_cell_iterator &cell) const;
+
+ // avoid -Woverloaded-virtual:
+ using Interface::properties_at_points;
+ };
+ }
+ }
+}
+
+#endif
diff --git a/source/particle/interpolator/bilinear_least_squares.cc b/source/particle/interpolator/bilinear_least_squares.cc
new file mode 100644
index 00000000000..b3e3c3955ed
--- /dev/null
+++ b/source/particle/interpolator/bilinear_least_squares.cc
@@ -0,0 +1,180 @@
+/*
+ Copyright (C) 2017 by the authors of the ASPECT code.
+
+ This file is part of ASPECT.
+
+ ASPECT is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2, or (at your option)
+ any later version.
+
+ ASPECT is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with ASPECT; see the file doc/COPYING. If not see
+ .
+ */
+
+#include
+
+#include
+#include
+#include
+
+namespace aspect
+{
+ namespace Particle
+ {
+ namespace Interpolator
+ {
+ template
+ std::vector >
+ BilinearLeastSquares::properties_at_points(const std::multimap > &particles,
+ const std::vector > &positions,
+ const ComponentMask &selected_properties,
+ const typename parallel::distributed::Triangulation::active_cell_iterator &cell) const
+ {
+ const unsigned int n_particle_properties = particles.begin()->second.get_properties().size();
+
+ const unsigned int property_index = selected_properties.first_selected_component(selected_properties.size());
+
+ AssertThrow(property_index != numbers::invalid_unsigned_int,
+ ExcMessage("Internal error: the particle property interpolator was "
+ "called without a specified component to interpolate."));
+
+ AssertThrow(dim == 2,
+ ExcMessage("Currently, the particle interpolator 'bilinear' is only supported for 2D models."));
+
+ AssertThrow(selected_properties.n_selected_components(n_particle_properties) == 1,
+ ExcNotImplemented("Interpolation of multiple components is not supported."));
+
+ const Point approximated_cell_midpoint = std::accumulate (positions.begin(), positions.end(), Point())
+ / static_cast (positions.size());
+
+ typename parallel::distributed::Triangulation::active_cell_iterator found_cell;
+
+ if (cell == typename parallel::distributed::Triangulation::active_cell_iterator())
+ {
+ // We can not simply use one of the points as input for find_active_cell_around_point
+ // because for vertices of mesh cells we might end up getting ghost_cells as return value
+ // instead of the local active cell. So make sure we are well in the inside of a cell.
+ Assert(positions.size() > 0,
+ ExcMessage("The particle property interpolator was not given any "
+ "positions to evaluate the particle cell_properties at."));
+
+
+ found_cell =
+ (GridTools::find_active_cell_around_point<> (this->get_mapping(),
+ this->get_triangulation(),
+ approximated_cell_midpoint)).first;
+ }
+ else
+ found_cell = cell;
+
+ const types::LevelInd cell_index = std::make_pair (found_cell->level(),found_cell->index());
+ const std::pair >::const_iterator,
+ typename std::multimap >::const_iterator> particle_range = particles.equal_range(cell_index);
+
+
+ std::vector > cell_properties(positions.size(),
+ std::vector(n_particle_properties,
+ numbers::signaling_nan()));
+
+ const unsigned int n_particles = std::distance(particle_range.first,particle_range.second);
+
+ AssertThrow(n_particles != 0,
+ ExcMessage("At least one cell contained no particles. The 'bilinear'"
+ "interpolation scheme does not support this case. "));
+
+
+ // Noticed that the size of matrix A is n_particles x matrix_dimension
+ // which usually is not a square matrix. Therefore, we solve Ax=r by
+ // solving A^TAx= A^Tr.
+ const unsigned int matrix_dimension = 4;
+ dealii::LAPACKFullMatrix A(n_particles, matrix_dimension);
+ Vector r(n_particles);
+ r = 0;
+
+ // Keep track of the local min and max in order to correct
+ // overshoot and undershoot of the interpolated particle property.
+ double max_value_for_particle_property = (particle_range.first)->second.get_properties()[property_index];
+ double min_value_for_particle_property = (particle_range.first)->second.get_properties()[property_index];
+
+ unsigned int index = 0;
+ const double cell_diameter = found_cell->diameter();
+ for (typename std::multimap >::const_iterator particle = particle_range.first;
+ particle != particle_range.second; ++particle, ++index)
+ {
+ const double particle_property_value = particle->second.get_properties()[property_index];
+ r[index] = particle_property_value;
+
+ const Point position = particle->second.get_location();
+ A(index,0) = 1;
+ A(index,1) = (position[0] - approximated_cell_midpoint[0])/cell_diameter;
+ A(index,2) = (position[1] - approximated_cell_midpoint[1])/cell_diameter;
+ A(index,3) = (position[0] - approximated_cell_midpoint[0]) * (position[1] - approximated_cell_midpoint[1])/std::pow(cell_diameter,2);
+
+ if (min_value_for_particle_property > particle_property_value)
+ min_value_for_particle_property = particle_property_value;
+ if (max_value_for_particle_property < particle_property_value)
+ max_value_for_particle_property = particle_property_value;
+ }
+
+ dealii::LAPACKFullMatrix B(matrix_dimension, matrix_dimension);
+ dealii::LAPACKFullMatrix B_inverse(B);
+
+ Vector c_ATr(matrix_dimension);
+ Vector c(matrix_dimension);
+
+ const double threshold = 1e-15;
+ unsigned int index_positions = 0;
+
+ // Matrix A can be rank deficient if it does not have full rank, therefore singular.
+ // To circumvent this issue, we solve A^TAx=A^Tr by using singular value
+ // decomposition (SVD).
+ A.Tmmult(B, A, false);
+ A.Tvmult(c_ATr,r);
+ B_inverse.compute_inverse_svd(threshold);
+ B_inverse.vmult(c, c_ATr);
+
+ for (typename std::vector >::const_iterator itr = positions.begin(); itr != positions.end(); ++itr, ++index_positions)
+ {
+ const Point support_point = *itr;
+ double interpolated_value = c[0] +
+ c[1]*(support_point[0] - approximated_cell_midpoint[0])/cell_diameter +
+ c[2]*(support_point[1] - approximated_cell_midpoint[1])/cell_diameter +
+ c[3]*(support_point[0] - approximated_cell_midpoint[0])*(support_point[1] - approximated_cell_midpoint[1])/std::pow(cell_diameter,2);
+
+ cell_properties[index_positions][property_index] = interpolated_value;
+
+ // Overshoot and undershoot correction of interpolated particle property.
+ if (interpolated_value > max_value_for_particle_property)
+ interpolated_value = max_value_for_particle_property;
+ else if (interpolated_value < min_value_for_particle_property)
+ interpolated_value = min_value_for_particle_property;
+ }
+ return cell_properties;
+ }
+ }
+ }
+}
+
+
+// explicit instantiations
+namespace aspect
+{
+ namespace Particle
+ {
+ namespace Interpolator
+ {
+ ASPECT_REGISTER_PARTICLE_INTERPOLATOR(BilinearLeastSquares,
+ "bilinear least squares",
+ "Interpolates particle properties onto a vector of points using a "
+ "bilinear least squares method. Currently only 2D models are "
+ "supported. Note that deal.II must be configured with BLAS/LAPACK.")
+ }
+ }
+}
diff --git a/tests/particle_interpolator_bilinear_solkz.cc b/tests/particle_interpolator_bilinear_solkz.cc
new file mode 100644
index 00000000000..3dd35c06451
--- /dev/null
+++ b/tests/particle_interpolator_bilinear_solkz.cc
@@ -0,0 +1,2 @@
+#include "../benchmarks/solkz/compositional_fields/solkz_compositional_fields.cc"
+
diff --git a/tests/particle_interpolator_bilinear_solkz.prm b/tests/particle_interpolator_bilinear_solkz.prm
new file mode 100644
index 00000000000..3e62de9aba2
--- /dev/null
+++ b/tests/particle_interpolator_bilinear_solkz.prm
@@ -0,0 +1,116 @@
+# A test based on the SolKz benchmark that performs bilinear
+# least squares interpolation from particle properties to
+# compositional fields that are flagged as 'particle' advected fields.
+
+set Additional shared libraries = ./libparticle_interpolator_bilinear_solkz.so
+
+############### Global parameters
+
+set Dimension = 2
+
+set Start time = 0
+set End time = 0
+
+set Output directory = output
+
+set Pressure normalization = volume
+
+############### Parameters describing the model
+
+subsection Geometry model
+ set Model name = box
+
+ subsection Box
+ set X extent = 1
+ set Y extent = 1
+ end
+end
+
+subsection Model settings
+ set Prescribed velocity boundary indicators =
+ set Tangential velocity boundary indicators = left, right, bottom, top
+ set Zero velocity boundary indicators =
+end
+
+subsection Material model
+ set Model name = SolKzCompositionalMaterial
+end
+
+subsection Gravity model
+ set Model name = vertical
+end
+
+############### Parameters describing the temperature field
+
+subsection Boundary temperature model
+ set Model name = box
+end
+
+subsection Initial temperature model
+ set Model name = perturbed box
+end
+
+############### Parameters describing the discretization
+
+subsection Discretization
+ set Stokes velocity polynomial degree = 2
+ set Use locally conservative discretization = false
+ set Use discontinuous composition discretization = true
+end
+
+subsection Mesh refinement
+ set Initial adaptive refinement = 0
+ set Initial global refinement = 5
+end
+
+############### Parameters describing the compositional field
+# Note: The compositional field is what drives the flow
+# in this example
+
+subsection Compositional fields
+ set Number of fields = 2
+ set Names of fields = density_comp, viscosity_comp
+ set Compositional field methods = particles, particles
+ set Mapped particle properties = density_comp:function[0], viscosity_comp:function[1]
+end
+
+subsection Initial composition model
+ set Model name = function
+ subsection Function
+ set Variable names = x,z
+ set Function constants = pi=3.1415926, eta_b=1e6
+ set Function expression = -1 * sin(2*z) * cos(3*pi*x); exp(log(eta_b)*z)
+ end
+end
+
+############### Parameters describing what to do with the solution
+
+subsection Postprocess
+ set List of postprocessors = particles, SolKzPostprocessor
+
+ subsection Particles
+ set Number of particles = 16384
+ set Time between data output = 0
+ set Data output format = vtu
+ set List of particle properties = function
+ set Integration scheme = rk2
+ set Interpolation scheme = cell average
+ set Maximum particles per cell = 16384
+ set Update ghost particles = true
+
+ subsection Function
+ set Number of components = 2
+ set Variable names = x, z
+ set Function constants = pi=3.1415926, eta_b=1e6
+ set Function expression = -1 * sin(2*z) * cos(3*pi*x); exp(log(eta_b)*z)
+ end
+
+ set Particle generator name = reference cell
+
+ subsection Generator
+ subsection Reference cell
+ set Number of particles per cell per direction = 4
+ end
+ end
+ end
+end
\ No newline at end of file
diff --git a/tests/particle_interpolator_bilinear_solkz/screen-output b/tests/particle_interpolator_bilinear_solkz/screen-output
new file mode 100644
index 00000000000..c69c4642fbf
--- /dev/null
+++ b/tests/particle_interpolator_bilinear_solkz/screen-output
@@ -0,0 +1,24 @@
+-----------------------------------------------------------------------------
+-----------------------------------------------------------------------------
+
+Loading shared library <./libparticle_interpolator_bilinear_solkz.so>
+
+Number of active cells: 1,024 (on 6 levels)
+Number of degrees of freedom: 32,196 (8,450+1,089+4,225+9,216+9,216)
+
+*** Timestep 0: t=0 years
+ Skipping temperature solve because RHS is zero.
+ Rebuilding Stokes preconditioner...
+ Solving Stokes system... 46+0 iterations.
+
+ Postprocessing:
+ Writing particle output: output-particle_interpolator_bilinear_solkz/particles/particles-00000
+ Errors u_L1, p_L1, u_L2, p_L2: 2.396693e-07, 3.813761e-04, 3.798509e-07, 1.430329e-03
+
+Termination requested by criterion: end time
+
+
++---------------------------------------------+------------+------------+
++---------------------------------+-----------+------------+------------+
++---------------------------------+-----------+------------+------------+
+
diff --git a/tests/particle_interpolator_bilinear_solkz/statistics b/tests/particle_interpolator_bilinear_solkz/statistics
new file mode 100644
index 00000000000..7c5077eba42
--- /dev/null
+++ b/tests/particle_interpolator_bilinear_solkz/statistics
@@ -0,0 +1,14 @@
+# 1: Time step number
+# 2: Time (years)
+# 3: Time step size (years)
+# 4: Number of mesh cells
+# 5: Number of Stokes degrees of freedom
+# 6: Number of temperature degrees of freedom
+# 7: Number of degrees of freedom for all compositions
+# 8: Iterations for temperature solver
+# 9: Iterations for Stokes solver
+# 10: Velocity iterations in Stokes preconditioner
+# 11: Schur complement iterations in Stokes preconditioner
+# 12: Number of advected particles
+# 13: Particle file name
+0 0.000000000000e+00 0.000000000000e+00 1024 9539 4225 18432 0 46 47 46 16384 output-particle_interpolator_bilinear_solkz/output-particle_interpolator_bilinear_solkz/particles/particles-00000