Qualify std functions everywhere

benchmarks/advection_in_annulus/advection_in_annulus.cc

@@ -64,10 +64,10 @@ namespace aspect
  const double theta = std::atan2(y,x);
  const double f_r = A*r + B/r;
  const double g_r = A*r/2 + B*std::log(r)/r + C/r;
- const double v_r = g_r*k*sin(k*theta);
- const double v_theta = f_r*cos(k*theta);
- const double v_x = cos(theta)*v_r - sin(theta)*v_theta;
- const double v_y = sin(theta)*v_r + cos(theta)*v_theta;
+ const double v_r = g_r*k*std::sin(k*theta);
+ const double v_theta = f_r*std::cos(k*theta);
+ const double v_x = std::cos(theta)*v_r - std::sin(theta)*v_theta;
+ const double v_y = std::sin(theta)*v_r + std::cos(theta)*v_theta;
  return Point<2> (v_x,v_y);
  }
 
@@ -82,7 +82,7 @@ namespace aspect
  const double f_r = 2*r + B/r;
  const double g_r = A*r/2 + B*std::log(r)/r + C/r;
  const double h_r=(2*g_r-f_r)/r;
- return k*h_r*sin(k*theta)+rho_0*gravity*(outer_radius-r);
+ return k*h_r*std::sin(k*theta)+rho_0*gravity*(outer_radius-r);
  }
 
  template <int dim>

benchmarks/annulus/plugin/annulus.cc

@@ -219,7 +219,7 @@ namespace aspect
  const double r = spherical_position[0];
  const double theta = spherical_position[1];
 
- const double forcing_term = m(r,k)*k*sin(k*(theta-phase(t))) + rho_0;
+ const double forcing_term = m(r,k)*k*std::sin(k*(theta-phase(t))) + rho_0;
  gravity_magnitude = forcing_term / Annulus_density(pos,k,t, transient);
  }
 

benchmarks/burstedde/burstedde.cc

@@ -192,7 +192,7 @@ namespace aspect
  const double x=p[0];
  const double y=p[1];
  const double z=p[2];
- const double mu=exp(1. - beta * (x*(1.-x)+y*(1.-y) + z*(1.-z)));
+ const double mu=std::exp(1. - beta * (x*(1.-x)+y*(1.-y) + z*(1.-z)));
 
  out.viscosities[i] = mu;
  out.thermal_conductivities[i] = 0.0;
@@ -352,7 +352,7 @@ namespace aspect
  const double x=pos[0];
  const double y=pos[1];
  const double z=pos[2];
- const double mu=exp(1. - beta * (x*(1.-x)+y*(1.-y) + z*(1.-z)));
+ const double mu=std::exp(1. - beta * (x*(1.-x)+y*(1.-y) + z*(1.-z)));
 
  const double dmudx=-beta*(1.-2.*x)*mu;
  const double dmudy=-beta*(1.-2.*y)*mu;

benchmarks/crameri_et_al/case_1/crameri_benchmark_1.cc

@@ -46,7 +46,7 @@ namespace aspect
  /**
  * Generate a coarse mesh for the geometry described by this class.
  * Makes perturbs the top boundary of the box with a function
- * of the form z' = amplitude * cos(order * x )
+ * of the form z' = amplitude * std::cos(order * x )
  */
  void create_coarse_mesh (parallel::distributed::Triangulation<dim> &coarse_grid) const override;
 

benchmarks/hollow_sphere/hollow_sphere.cc

@@ -73,25 +73,25 @@ namespace aspect
 
  if (mmm == -1)
  {
- alpha=-gammma*(pow(R2,3)-pow(R1,3))/(pow(R2,3)*log(R1)-pow(R1,3)*log(R2));
- beta=-3*gammma*(log(R2)-log(R1))/(pow(R1,3)*log(R2)-pow(R2,3)*log(R1)) ;
+ alpha=-gammma*(std::pow(R2,3)-std::pow(R1,3))/(std::pow(R2,3)*std::log(R1)-std::pow(R1,3)*std::log(R2));
+ beta=-3*gammma*(std::log(R2)-std::log(R1))/(std::pow(R1,3)*std::log(R2)-std::pow(R2,3)*std::log(R1)) ;
  fr=alpha/(r*r)+beta*r;
- gr=-2/(r*r)*(alpha*log(r)+beta/3*pow(r,3)+gammma);
+ gr=-2/(r*r)*(alpha*std::log(r)+beta/3*std::pow(r,3)+gammma);
  }
  else
  {
- alpha=gammma*(mmm+1)*(pow(R1,-3)-pow(R2,-3))/(pow(R1,-mmm-4)-pow(R2,-mmm-4));
- beta=-3*gammma*(pow(R1,mmm+1)-pow(R2,mmm+1))/(pow(R1,mmm+4)-pow(R2,mmm+4));
- fr=alpha/pow(r,mmm+3)+beta*r;
- gr=-2/(r*r)*(-alpha/(mmm+1)*pow(r,-mmm-1)+beta/3*pow(r,3)+gammma);
+ alpha=gammma*(mmm+1)*(std::pow(R1,-3)-std::pow(R2,-3))/(std::pow(R1,-mmm-4)-std::pow(R2,-mmm-4));
+ beta=-3*gammma*(std::pow(R1,mmm+1)-std::pow(R2,mmm+1))/(std::pow(R1,mmm+4)-std::pow(R2,mmm+4));
+ fr=alpha/std::pow(r,mmm+3)+beta*r;
+ gr=-2/(r*r)*(-alpha/(mmm+1)*std::pow(r,-mmm-1)+beta/3*std::pow(r,3)+gammma);
  }
 
- const double v_r =gr*cos(theta);
- const double v_theta=fr*sin(theta);
- const double v_phi =fr*sin(theta);
- const double v_x=sin(theta)*cos(phi)*v_r + cos(theta)*cos(phi)*v_theta-sin(phi)*v_phi;
- const double v_y=sin(theta)*sin(phi)*v_r + cos(theta)*sin(phi)*v_theta+cos(phi)*v_phi;
- const double v_z=cos(theta)*v_r - sin(theta)*v_theta;
+ const double v_r =gr*std::cos(theta);
+ const double v_theta=fr*std::sin(theta);
+ const double v_phi =fr*std::sin(theta);
+ const double v_x=std::sin(theta)*std::cos(phi)*v_r + std::cos(theta)*std::cos(phi)*v_theta-std::sin(phi)*v_phi;
+ const double v_y=std::sin(theta)*std::sin(phi)*v_r + std::cos(theta)*std::sin(phi)*v_theta+std::cos(phi)*v_phi;
+ const double v_z=std::cos(theta)*v_r - std::sin(theta)*v_theta;
 
  // create a Point<3> (because it has a constructor that takes
  // three doubles) and return it (it automatically converts to
@@ -115,21 +115,21 @@ namespace aspect
  if (mmm == -1)
  {
  mur=mu0;
- alpha=-gammma*(pow(R2,3)-pow(R1,3))/(pow(R2,3)*log(R1)-pow(R1,3)*log(R2));
- beta=-3*gammma*(log(R2)-log(R1))/(pow(R1,3)*log(R2)-pow(R2,3)*log(R1)) ;
- gr=-2/(r*r)*(alpha*log(r)+beta/3*pow(r,3)+gammma);
+ alpha=-gammma*(std::pow(R2,3)-std::pow(R1,3))/(std::pow(R2,3)*std::log(R1)-std::pow(R1,3)*std::log(R2));
+ beta=-3*gammma*(std::log(R2)-std::log(R1))/(std::pow(R1,3)*std::log(R2)-std::pow(R2,3)*std::log(R1)) ;
+ gr=-2/(r*r)*(alpha*std::log(r)+beta/3*std::pow(r,3)+gammma);
  hr=2/r*gr*mur;
  }
  else
  {
- mur=mu0*pow(r,mmm+1);
- alpha=gammma*(mmm+1)*(pow(R1,-3)-pow(R2,-3))/(pow(R1,-mmm-4)-pow(R2,-mmm-4));
- beta=-3*gammma*(pow(R1,mmm+1)-pow(R2,mmm+1))/(pow(R1,mmm+4)-pow(R2,mmm+4));
- gr=-2/(r*r)*(-alpha/(mmm+1)*pow(r,-mmm-1)+beta/3*pow(r,3)+gammma);
+ mur=mu0*std::pow(r,mmm+1);
+ alpha=gammma*(mmm+1)*(std::pow(R1,-3)-std::pow(R2,-3))/(std::pow(R1,-mmm-4)-std::pow(R2,-mmm-4));
+ beta=-3*gammma*(std::pow(R1,mmm+1)-std::pow(R2,mmm+1))/(std::pow(R1,mmm+4)-std::pow(R2,mmm+4));
+ gr=-2/(r*r)*(-alpha/(mmm+1)*std::pow(r,-mmm-1)+beta/3*std::pow(r,3)+gammma);
  hr=(mmm+3)/r*gr*mur;
  }
 
- return hr*cos(theta) + rho_0 * gravity * (R2 - r);
+ return hr*std::cos(theta) + rho_0 * gravity * (R2 - r);
  }
 
  template <int dim>
@@ -147,20 +147,20 @@ namespace aspect
 
  if (mmm == -1)
  {
- alpha=-gammma*(pow(R2,3)-pow(R1,3))/(pow(R2,3)*log(R1)-pow(R1,3)*log(R2));
- beta=-3*gammma*(log(R2)-log(R1))/(pow(R1,3)*log(R2)-pow(R2,3)*log(R1)) ;
+ alpha=-gammma*(std::pow(R2,3)-std::pow(R1,3))/(std::pow(R2,3)*std::log(R1)-std::pow(R1,3)*std::log(R2));
+ beta=-3*gammma*(std::log(R2)-std::log(R1))/(std::pow(R1,3)*std::log(R2)-std::pow(R2,3)*std::log(R1)) ;
  fr=alpha/(r*r)+beta*r;
- gr=-2/(r*r)*(alpha*log(r)+beta/3*pow(r,3)+gammma);
+ gr=-2/(r*r)*(alpha*std::log(r)+beta/3*std::pow(r,3)+gammma);
  }
  else
  {
- alpha=gammma*(mmm+1)*(pow(R1,-3)-pow(R2,-3))/(pow(R1,-mmm-4)-pow(R2,-mmm-4));
- beta=-3*gammma*(pow(R1,mmm+1)-pow(R2,mmm+1))/(pow(R1,mmm+4)-pow(R2,mmm+4));
- fr=alpha/pow(r,mmm+3)+beta*r;
- gr=-2/(r*r)*(-alpha/(mmm+1)*pow(r,-mmm-1)+beta/3*pow(r,3)+gammma);
+ alpha=gammma*(mmm+1)*(std::pow(R1,-3)-std::pow(R2,-3))/(std::pow(R1,-mmm-4)-std::pow(R2,-mmm-4));
+ beta=-3*gammma*(std::pow(R1,mmm+1)-std::pow(R2,mmm+1))/(std::pow(R1,mmm+4)-std::pow(R2,mmm+4));
+ fr=alpha/std::pow(r,mmm+3)+beta*r;
+ gr=-2/(r*r)*(-alpha/(mmm+1)*std::pow(r,-mmm-1)+beta/3*std::pow(r,3)+gammma);
  }
 
- return -(6.*gr + 4.*fr) * cos(theta) * mu0 / r;
+ return -(6.*gr + 4.*fr) * std::cos(theta) * mu0 / r;
  }
 
 
@@ -366,7 +366,7 @@ namespace aspect
  const Point<dim> &pos = in.position[i];
  const std::array<double,dim> spos = aspect::Utilities::Coordinates::cartesian_to_spherical_coordinates(pos);
  const double r = spos[0];
- const double mu = pow(r,mmm+1);
+ const double mu = std::pow(r,mmm+1);
  out.viscosities[i] = mu;
 
  const double theta=spos[2];
@@ -380,15 +380,15 @@ namespace aspect
 
  if (mmm == -1)
  {
- alpha = -gammma*(pow(R2,3)-pow(R1,3))/(pow(R2,3)*log(R1)-pow(R1,3)*log(R2));
- beta = -3*gammma*(log(R2)-log(R1))/(pow(R1,3)*log(R2)-pow(R2,3)*log(R1)) ;
- rho = -(alpha/pow(r,4)*(8*log(r)-6) + 8./3.*beta/r+8*gammma/pow(r,4))*cos(theta) + rho_0;
+ alpha = -gammma*(std::pow(R2,3)-std::pow(R1,3))/(std::pow(R2,3)*std::log(R1)-std::pow(R1,3)*std::log(R2));
+ beta = -3*gammma*(std::log(R2)-std::log(R1))/(std::pow(R1,3)*std::log(R2)-std::pow(R2,3)*std::log(R1)) ;
+ rho = -(alpha/std::pow(r,4)*(8*std::log(r)-6) + 8./3.*beta/r+8*gammma/std::pow(r,4))*std::cos(theta) + rho_0;
  }
  else
  {
- alpha=gammma*(mmm+1)*(pow(R1,-3)-pow(R2,-3))/(pow(R1,-mmm-4)-pow(R2,-mmm-4));
- beta=-3*gammma*(pow(R1,mmm+1)-pow(R2,mmm+1))/(pow(R1,mmm+4)-pow(R2,mmm+4));
- rho= -(2*alpha*pow(r,-4)*(mmm+3)/(mmm+1)*(mmm-1)-2*beta/3*(mmm-1)*(mmm+3)*pow(r,mmm)-mmm*(mmm+5)*2*gammma*pow(r,mmm-3) )*cos(theta) + rho_0;
+ alpha=gammma*(mmm+1)*(std::pow(R1,-3)-std::pow(R2,-3))/(std::pow(R1,-mmm-4)-std::pow(R2,-mmm-4));
+ beta=-3*gammma*(std::pow(R1,mmm+1)-std::pow(R2,mmm+1))/(std::pow(R1,mmm+4)-std::pow(R2,mmm+4));
+ rho= -(2*alpha*std::pow(r,-4)*(mmm+3)/(mmm+1)*(mmm-1)-2*beta/3*(mmm-1)*(mmm+3)*std::pow(r,mmm)-mmm*(mmm+5)*2*gammma*std::pow(r,mmm-3) )*std::cos(theta) + rho_0;
  }
 
  out.densities[i] = rho;

benchmarks/king2dcompressible/plugin/code.cc

@@ -72,13 +72,13 @@ namespace aspect
 
  const double depth = 1.0-position(dim-1);
 
- out.viscosities[i] = ((Di==0.0)?1.0:Di)/Ra*exp(-b*(temperature- this->get_adiabatic_conditions().temperature(position))+c*depth);
+ out.viscosities[i] = ((Di==0.0)?1.0:Di)/Ra*std::exp(-b*(temperature- this->get_adiabatic_conditions().temperature(position))+c*depth);
 
  out.specific_heat[i] = reference_specific_heat;
  out.thermal_conductivities[i] = 1.0;
  out.thermal_expansion_coefficients[i] = (Di==0.0)?1.0:Di;
 
- double rho = reference_rho * exp(depth * Di/gamma);
+ double rho = reference_rho * std::exp(depth * Di/gamma);
  rho *= 1.0 - out.thermal_expansion_coefficients[i] * (temperature - this->get_adiabatic_conditions().temperature(position))
  + (tala?0.0:1.0)*Di*gamma
  * (pressure - this->get_adiabatic_conditions().pressure(position));

benchmarks/layeredflow/layeredflow.cc

@@ -50,8 +50,8 @@ namespace aspect
  const double yplus=(1.+y0)/beta;
  const double yminus=(1.-y0)/beta;
  const double C1 = 2.*numbers::PI /
- ( beta*std::log(beta*beta+pow(1.+y0,2.0))-2.*(1.+y0)*std::atan(yplus)
- - beta*std::log(beta*beta+pow(1.-y0,2.0))+2.*(1.-y0)*std::atan(yminus)
+ ( beta*std::log(beta*beta+std::pow(1.+y0,2.0))-2.*(1.+y0)*std::atan(yplus)
+ - beta*std::log(beta*beta+std::pow(1.-y0,2.0))+2.*(1.-y0)*std::atan(yminus)
  + 2.*numbers::PI*(1.+2.*epsilon) );
 
  const double C2 = ( beta*std::log( beta*beta+Utilities::fixed_power<2>(1+y0) )- 2.*(1.+y0)*std::atan(yplus)

benchmarks/operator_splitting/advection_reaction/advection_reaction.cc

@@ -126,7 +126,7 @@ namespace aspect
  for (unsigned int q=0; q < in.n_evaluation_points(); ++q)
  {
  // dC/dt = - z * lambda * C
- const double decay_constant = half_life > 0.0 ? log(2.0) / half_life : 0.0;
+ const double decay_constant = half_life > 0.0 ? std::log(2.0) / half_life : 0.0;
  const double z = in.position[q](1);
  reaction_out->reaction_rates[q][0] = - decay_constant * z / time_scale * in.composition[q][0];
  }
@@ -225,7 +225,7 @@ namespace aspect
  for (unsigned int q=0; q<heating_model_outputs.heating_source_terms.size(); ++q)
  {
  // dC/dt = - z * lambda * C
- const double decay_constant = half_life > 0.0 ? log(2.0) / half_life : 0.0;
+ const double decay_constant = half_life > 0.0 ? std::log(2.0) / half_life : 0.0;
  const double z = in.position[q](1);
  heating_model_outputs.rates_of_temperature_change[q] = - decay_constant * z / time_scale * in.composition[q][0];
 
@@ -290,8 +290,8 @@ namespace aspect
  values[0] = 0.0; // velocity x
  values[1] = 0.0; // velocity z
  values[2] = 0.0; // pressure
- values[3] = sin(2*numbers::PI*(x-t*0.01))*exp(-log(2.0)/10.0*z*t); // temperature
- values[4] = sin(2*numbers::PI*(x-t*0.01))*exp(-log(2.0)/10.0*z*t); // composition
+ values[3] = std::sin(2*numbers::PI*(x-t*0.01))*std::exp(-std::log(2.0)/10.0*z*t); // temperature
+ values[4] = std::sin(2*numbers::PI*(x-t*0.01))*std::exp(-std::log(2.0)/10.0*z*t); // composition
  }
  };
 
@@ -324,7 +324,7 @@ namespace aspect
 
  double t = this->get_time() / time_scale;
 
- return sin(2*numbers::PI*(x-t*0.01))*exp(-log(2.0)/10.0*z*t);
+ return std::sin(2*numbers::PI*(x-t*0.01))*std::exp(-std::log(2.0)/10.0*z*t);
  }
 
 
@@ -358,7 +358,7 @@ namespace aspect
 
  double t = this->get_time() / time_scale;
 
- return sin(2*numbers::PI*(x-t*0.01))*exp(-log(2.0)/10.0*z*t);
+ return std::sin(2*numbers::PI*(x-t*0.01))*std::exp(-std::log(2.0)/10.0*z*t);
  }
 
  template <int dim>

benchmarks/operator_splitting/exponential_decay/exponential_decay.cc

@@ -126,7 +126,7 @@ namespace aspect
  for (unsigned int q=0; q < in.n_evaluation_points(); ++q)
  {
  // dC/dt = - lambda * C
- const double decay_constant = half_life > 0.0 ? log(2.0) / half_life : 0.0;
+ const double decay_constant = half_life > 0.0 ? std::log(2.0) / half_life : 0.0;
  reaction_out->reaction_rates[q][0] = - decay_constant / time_scale * in.composition[q][0];
  }
  }
@@ -225,7 +225,7 @@ namespace aspect
  for (unsigned int q=0; q<heating_model_outputs.heating_source_terms.size(); ++q)
  {
  // dC/dt = - lambda * C
- const double decay_constant = half_life > 0.0 ? log(2.0) / half_life : 0.0;
+ const double decay_constant = half_life > 0.0 ? std::log(2.0) / half_life : 0.0;
  heating_model_outputs.rates_of_temperature_change[q] = - decay_constant / time_scale * in.composition[q][0];
 
  heating_model_outputs.heating_source_terms[q] = 0.0;
@@ -285,8 +285,8 @@ namespace aspect
  values[0] = 0.0; // velocity x
  values[1] = 0.0; // velocity z
  values[2] = 0.0; // pressure
- values[3] = exp(-log(2.0)/10.0*this->get_time()); // temperature
- values[4] = exp(-log(2.0)/10.0*this->get_time()); // composition
+ values[3] = std::exp(-std::log(2.0)/10.0*this->get_time()); // temperature
+ values[4] = std::exp(-std::log(2.0)/10.0*this->get_time()); // composition
  }
  };
 

benchmarks/shear_bands/shear_bands.cc

@@ -86,7 +86,7 @@ namespace aspect
 
  double reference_darcy_coefficient () const override
  {
- return reference_permeability * pow(0.01, permeability_exponent) / eta_f;
+ return reference_permeability * std::pow(0.01, permeability_exponent) / eta_f;
  }
 
 
@@ -150,7 +150,7 @@ namespace aspect
  {
  double porosity = std::max(in.composition[i][porosity_idx],1e-4);
 
- melt_out->compaction_viscosities[i] = xi_0 * pow(porosity/background_porosity,-compaction_viscosity_exponent);
+ melt_out->compaction_viscosities[i] = xi_0 * std::pow(porosity/background_porosity,-compaction_viscosity_exponent);
  melt_out->fluid_viscosities[i]= eta_f;
  melt_out->permeabilities[i]= reference_permeability * std::pow(porosity,permeability_exponent);
  melt_out->fluid_densities[i]= reference_rho_f;
@@ -534,8 +534,8 @@ namespace aspect
  PlaneWaveMeltBandsInitialCondition<dim>::
  initial_composition (const Point<dim> &position, const unsigned int /*n_comp*/) const
  {
- return background_porosity * (1.0 + amplitude * cos(wave_number*position[0]*sin(initial_band_angle)
- + wave_number*position[1]*cos(initial_band_angle)));
+ return background_porosity * (1.0 + amplitude * std::cos(wave_number*position[0]*std::sin(initial_band_angle)
+  + wave_number*position[1]*std::cos(initial_band_angle)));
  }
 
 
@@ -736,7 +736,7 @@ namespace aspect
  const double min_velocity = bm.boundary_velocity(lower_boundary, lower_boundary_point).norm();
 
  const double strain_rate = 0.5 * (max_velocity + min_velocity) / this->get_geometry_model().maximal_depth();
- const double theta = std::atan(std::sin(initial_band_angle) / (std::cos(initial_band_angle) - time * strain_rate/sqrt(2.0) * std::sin(initial_band_angle)));
+ const double theta = std::atan(std::sin(initial_band_angle) / (std::cos(initial_band_angle) - time * strain_rate/numbers::SQRT2 * std::sin(initial_band_angle)));
  const double analytical_growth_rate = - eta_0 / (xi_0 + 4.0/3.0 * eta_0) * alpha * (1.0 - background_porosity)
  * 2.0 * strain_rate * std::sin(2.0 * theta);
 

benchmarks/solitary_wave/solitary_wave.cc

@@ -393,7 +393,7 @@ namespace aspect
  initial_composition_manager->template
  get_matching_initial_composition_model<SolitaryWaveInitialCondition<dim>>();
 
- return reference_permeability * pow(initial_composition.get_background_porosity(), 3.0) / eta_f;
+ return reference_permeability * std::pow(initial_composition.get_background_porosity(), 3.0) / eta_f;
 
  }
 

benchmarks/solubility/plugin/solubility.cc

@@ -156,7 +156,7 @@ namespace aspect
  for (unsigned int q=0; q<out.n_evaluation_points(); ++q)
  {
  const double porosity = std::max(in.composition[q][porosity_idx],0.0);
- out.viscosities[q] *= (1.0 - porosity) * exp(- alpha_phi * porosity);
+ out.viscosities[q] *= (1.0 - porosity) * std::exp(- alpha_phi * porosity);
  }
  }
 

benchmarks/tangurnis/code/tangurnis.cc

@@ -90,10 +90,10 @@ namespace aspect
 
  const Point<dim> &pos = in.position[i];
  const double depth = 1.0 - pos[dim-1];
- const double temperature = sin(numbers::PI*pos(dim-1))*cos(numbers::PI*wavenumber*pos(0));
+ const double temperature = std::sin(numbers::PI*pos(dim-1))*std::cos(numbers::PI*wavenumber*pos(0));
 
- out.viscosities[i] = ( Di==0.0 ? 1.0 : Di ) * exp( a * depth );
- out.densities[i] = ( Di==0.0 ? 1.0 : Di ) * (-1.0 * temperature ) * exp( Di/gamma * (depth) );
+ out.viscosities[i] = ( Di==0.0 ? 1.0 : Di ) * std::exp( a * depth );
+ out.densities[i] = ( Di==0.0 ? 1.0 : Di ) * (-1.0 * temperature ) * std::exp( Di/gamma * (depth) );
  out.specific_heat[i] = 1.0;
  out.thermal_conductivities[i] = 1.0;
  out.thermal_expansion_coefficients[i] = ( Di==0.0 ) ? 1.0 : Di;
@@ -293,7 +293,7 @@ namespace aspect
  const Point<dim> &position) const override
  {
  double wavenumber=1;
- return sin(numbers::PI*position(dim-1))*cos(numbers::PI*wavenumber*position(0));
+ return std::sin(numbers::PI*position(dim-1))*std::cos(numbers::PI*wavenumber*position(0));
  }
 
  double minimal_temperature (const std::set<types::boundary_id> &fixed_boundary_ids) const override;