ASPECT
interface.h
Go to the documentation of this file.
1 /*
2  Copyright (C) 2011 - 2020 by the authors of the ASPECT code.
3 
4  This file is part of ASPECT.
5 
6  ASPECT is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation; either version 2, or (at your option)
9  any later version.
10 
11  ASPECT is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with ASPECT; see the file LICENSE. If not see
18  <http://www.gnu.org/licenses/>.
19 */
20 
21 #ifndef _aspect_material_model_interface_h
22 #define _aspect_material_model_interface_h
23 
24 #include <aspect/global.h>
25 #include <aspect/plugins.h>
27 
28 #include <deal.II/base/point.h>
29 
30 // Work around an incorrect instantiation in qprojector.h of deal.II 9.2.0,
31 // which requires including qprojector.h before quadrature.h (and not
32 // after). This file doesn't actually need qprojector.h, so the include can be
33 // removed when we require 9.3.. For more info see
34 // https://github.com/geodynamics/aspect/issues/3728
35 #include <deal.II/base/quadrature.h>
36 
37 #include <deal.II/base/symmetric_tensor.h>
38 #include <deal.II/base/parameter_handler.h>
39 #include <deal.II/dofs/dof_handler.h>
40 #include <deal.II/dofs/dof_accessor.h>
41 #include <deal.II/fe/mapping.h>
42 #include <deal.II/fe/fe_values.h>
43 #include <deal.II/fe/component_mask.h>
44 #include <deal.II/numerics/data_postprocessor.h>
45 #include <deal.II/base/signaling_nan.h>
46 
47 namespace aspect
48 {
49  template <int dim>
50  struct Introspection;
51 
52 
60  namespace MaterialModel
61  {
62  using namespace dealii;
63 
68  namespace NonlinearDependence
69  {
92  {
94 
95  none = 1,
97  pressure = 4,
100 
102  };
103 
104 
109  const Dependence d2)
110  {
111  return Dependence(static_cast<int>(d1) | static_cast<int>(d2));
112  }
113 
115  const Dependence d2)
116  {
117  d1 = (d1 | d2);
118  return d1;
119  }
120 
126  {
132 
138 
144 
150 
156 
162  ModelDependence ();
163  };
164 
173  bool
174  identifies_single_variable(const Dependence dependence);
175  }
176 
181  namespace MaterialProperties
182  {
192  enum Property
193  {
195 
196  none = 1,
198  density = 4,
206 
209  specific_heat |
214  viscosity |
217  };
218 
223  inline Property operator | (const Property d1,
224  const Property d2)
225  {
226  return Property(static_cast<int>(d1) | static_cast<int>(d2));
227  }
228 
230  const Property d2)
231  {
232  return (d1 | d2);
233  }
234  }
235 
236  template <int dim> class AdditionalMaterialInputs;
237 
244  template <int dim>
246  {
257  MaterialModelInputs(const unsigned int n_points,
258  const unsigned int n_comp);
259 
275  MaterialModelInputs(const DataPostprocessorInputs::Vector<dim> &input_data,
276  const Introspection<dim> &introspection,
277  const bool compute_strain_rate = true);
278 
279 
299  MaterialModelInputs(const FEValuesBase<dim,dim> &fe_values,
300  const typename DoFHandler<dim>::active_cell_iterator &cell,
301  const Introspection<dim> &introspection,
302  const LinearAlgebra::BlockVector &solution_vector,
303  const bool compute_strain_rate = true);
304 
319 
324 
329  MaterialModelInputs &operator= (const MaterialModelInputs &source) = delete;
330 
334  MaterialModelInputs &operator= (MaterialModelInputs &&) = default;
335 
341  void reinit(const FEValuesBase<dim,dim> &fe_values,
342  const typename DoFHandler<dim>::active_cell_iterator &cell,
343  const Introspection<dim> &introspection,
344  const LinearAlgebra::BlockVector &solution_vector,
345  const bool compute_strain_rate = true);
346 
351  unsigned int n_evaluation_points() const;
352 
359  bool requests_property(const MaterialProperties::Property &property) const;
360 
365  std::vector<Point<dim>> position;
366 
370  std::vector<double> temperature;
371 
375  std::vector<double> pressure;
376 
381  std::vector<Tensor<1,dim>> pressure_gradient;
382 
390  std::vector<Tensor<1,dim>> velocity;
391 
397  std::vector<std::vector<double>> composition;
398 
411  std::vector<SymmetricTensor<2,dim>> strain_rate;
412 
428  typename DoFHandler<dim>::active_cell_iterator current_cell;
429 
438 
444  std::vector<std::unique_ptr<AdditionalMaterialInputs<dim>>> additional_inputs;
445 
452  template <class AdditionalInputType>
453  AdditionalInputType *get_additional_input();
454 
458  template <class AdditionalInputType>
459  const AdditionalInputType *get_additional_input() const;
460  };
461 
462 
463  template <int dim> class AdditionalMaterialOutputs;
464 
465 
472  template <int dim>
474  {
485  MaterialModelOutputs (const unsigned int n_points,
486  const unsigned int n_comp);
487 
502 
507 
512  MaterialModelOutputs &operator= (const MaterialModelOutputs &source) = delete;
513 
517  MaterialModelOutputs &operator= (MaterialModelOutputs &&) = default;
518 
523  unsigned int n_evaluation_points() const;
524 
528  std::vector<double> viscosities;
529 
533  std::vector<double> densities;
534 
539  std::vector<double> thermal_expansion_coefficients;
540 
544  std::vector<double> specific_heat;
545 
549  std::vector<double> thermal_conductivities;
550 
555  std::vector<double> compressibilities;
556 
562  std::vector<double> entropy_derivative_pressure;
563 
570  std::vector<double> entropy_derivative_temperature;
571 
608  std::vector<std::vector<double>> reaction_terms;
609 
615  std::vector<std::unique_ptr<AdditionalMaterialOutputs<dim>>> additional_outputs;
616 
624  template <class AdditionalOutputType>
625  AdditionalOutputType *get_additional_output();
626 
630  template <class AdditionalOutputType>
631  const AdditionalOutputType *get_additional_output() const;
632 
637  void move_additional_outputs_from(MaterialModelOutputs<dim> &other);
638  };
639 
640 
641 
656  namespace MaterialAveraging
657  {
697  {
707  };
708 
709 
716  std::string get_averaging_operation_names ();
717 
723  AveragingOperation parse_averaging_operation_name (const std::string &s);
724 
730  template <int dim>
731  void average (const AveragingOperation operation,
732  const typename DoFHandler<dim>::active_cell_iterator &cell,
733  const Quadrature<dim> &quadrature_formula,
734  const Mapping<dim> &mapping,
735  MaterialModelOutputs<dim> &values_out);
736 
742  void average_property (const AveragingOperation operation,
743  const FullMatrix<double> &projection_matrix,
744  const FullMatrix<double> &expansion_matrix,
745  std::vector<double> &values_out);
746  }
747 
748 
767  template <int dim>
769  {
770  public:
771  virtual ~AdditionalMaterialInputs() = default;
772 
777  virtual void
778  fill (const LinearAlgebra::BlockVector &solution,
779  const FEValuesBase<dim> &fe_values,
780  const Introspection<dim> &introspection) = 0;
781  };
782 
783 
810  template <int dim>
812  {
813  public:
818  virtual ~AdditionalMaterialOutputs() = default;
819 
820  virtual void average (const MaterialAveraging::AveragingOperation /*operation*/,
821  const FullMatrix<double> &/*projection_matrix*/,
822  const FullMatrix<double> &/*expansion_matrix*/)
823  {}
824  };
825 
826 
848  template <int dim>
850  {
851  public:
860  NamedAdditionalMaterialOutputs(const std::vector<std::string> &output_names);
861 
872  NamedAdditionalMaterialOutputs(const std::vector<std::string> &output_names,
873  const unsigned int n_points);
874 
878  ~NamedAdditionalMaterialOutputs() override;
879 
884  const std::vector<std::string> &get_names() const;
885 
890  virtual std::vector<double> get_nth_output(const unsigned int idx) const;
891 
893  const FullMatrix<double> &/*projection_matrix*/,
894  const FullMatrix<double> &/*expansion_matrix*/) override
895  {}
896 
897 
902  std::vector<std::vector<double>> output_values;
903 
904  private:
905  const std::vector<std::string> names;
906  };
907 
908 
914  template <int dim>
916  {
917  public:
918  SeismicAdditionalOutputs(const unsigned int n_points);
919 
920  std::vector<double> get_nth_output(const unsigned int idx) const override;
921 
927  std::vector<double> vs;
928 
934  std::vector<double> vp;
935  };
936 
937 
958  template <int dim>
960  {
961  public:
962  ReactionRateOutputs (const unsigned int n_points,
963  const unsigned int n_comp);
964 
965  std::vector<double> get_nth_output(const unsigned int idx) const override;
966 
974  std::vector<std::vector<double>> reaction_rates;
975  };
976 
977 
978 
984  template <int dim>
986  {
987  public:
988  PhaseOutputs(const unsigned int n_points);
989  };
990 
991 
992 
1011  template <int dim>
1013  {
1014  public:
1015  PrescribedFieldOutputs (const unsigned int n_points,
1016  const unsigned int n_comp);
1017 
1018  std::vector<double> get_nth_output(const unsigned int idx) const override;
1019 
1027  std::vector<std::vector<double>> prescribed_field_outputs;
1028  };
1029 
1043  template <int dim>
1045  {
1046  public:
1047  PrescribedTemperatureOutputs (const unsigned int n_points);
1048 
1049  std::vector<double> get_nth_output(const unsigned int idx) const override;
1050 
1058  std::vector<double> prescribed_temperature_outputs;
1059  };
1060 
1067  template <int dim>
1069  {
1070  public:
1071  AdditionalMaterialOutputsStokesRHS(const unsigned int n_points)
1072  : rhs_u(n_points), rhs_p(n_points), rhs_melt_pc(n_points)
1073  {}
1074 
1076  {}
1077 
1079  const FullMatrix<double> &/*projection_matrix*/,
1080  const FullMatrix<double> &/*expansion_matrix*/) override
1081  {
1082  // TODO: not implemented
1083  }
1084 
1090  std::vector<Tensor<1,dim>> rhs_u;
1091 
1096  std::vector<double> rhs_p;
1097 
1102  std::vector<double> rhs_melt_pc;
1103  };
1104 
1105 
1106 
1122  template <int dim>
1124  {
1125  public:
1129  explicit PrescribedPlasticDilation (const unsigned int n_points);
1130 
1134  virtual std::vector<double> get_nth_output(const unsigned int idx) const;
1135 
1140  std::vector<double> dilation;
1141  };
1142 
1143 
1144 
1151  template <int dim>
1153  {
1154  public:
1155  ElasticOutputs(const unsigned int n_points)
1156  : elastic_force(n_points, numbers::signaling_nan<Tensor<2,dim>>() )
1157  {}
1158 
1159  ~ElasticOutputs() override
1160  {}
1161 
1163  const FullMatrix<double> &/*projection_matrix*/,
1164  const FullMatrix<double> &/*expansion_matrix*/) override
1165  {
1166  AssertThrow(operation == MaterialAveraging::AveragingOperation::none,ExcNotImplemented());
1167  return;
1168  }
1169 
1175  std::vector<Tensor<2,dim>> elastic_force;
1176  };
1177 
1178 
1179 
1195  template <int dim>
1196  class Interface
1197  {
1198  public:
1213 
1218  virtual ~Interface() = default;
1219 
1225  virtual
1226  void
1227  initialize ();
1228 
1233  virtual void update ();
1234 
1247  get_model_dependence () const;
1248 
1258  virtual bool is_compressible () const = 0;
1285  virtual double reference_viscosity () const = 0;
1295  virtual
1296  void evaluate (const MaterialModel::MaterialModelInputs<dim> &in,
1308  static
1309  void
1310  declare_parameters (ParameterHandler &prm);
1311 
1318  virtual
1319  void
1320  parse_parameters (ParameterHandler &prm);
1330  virtual
1331  void
1332  create_additional_named_outputs (MaterialModelOutputs &outputs) const;
1333 
1334 
1342  virtual
1343  void
1344  fill_additional_material_model_inputs(MaterialModel::MaterialModelInputs<dim> &input,
1345  const LinearAlgebra::BlockVector &solution,
1346  const FEValuesBase<dim> &fe_values,
1347  const Introspection<dim> &introspection) const;
1348 
1349  protected:
1366  };
1367 
1383  template <int dim>
1384  void
1385  register_material_model (const std::string &name,
1386  const std::string &description,
1387  void (*declare_parameters_function) (ParameterHandler &),
1388  Interface<dim> *(*factory_function) ());
1389 
1400  template <int dim>
1401  Interface<dim> *
1402  create_material_model (const std::string &model_name);
1403 
1404 
1416  template <int dim>
1417  Interface<dim> *
1418  create_material_model (ParameterHandler &prm);
1419 
1420 
1427  template <int dim>
1428  std::string
1430 
1431 
1437  template <int dim>
1438  void
1439  declare_parameters (ParameterHandler &prm);
1440 
1441 
1442 
1452  template <int dim>
1453  void
1454  write_plugin_graph (std::ostream &output_stream);
1455 
1456 
1457 
1458 // --------------------- template function definitions ----------------------------------
1459 
1460  template <int dim>
1461  template <class AdditionalInputType>
1463  {
1464  for (unsigned int i=0; i<additional_inputs.size(); ++i)
1465  {
1466  AdditionalInputType *result = dynamic_cast<AdditionalInputType *> (additional_inputs[i].get());
1467  if (result)
1468  return result;
1469  }
1470  return nullptr;
1471  }
1472 
1473 
1474  template <int dim>
1475  template <class AdditionalInputType>
1476  const AdditionalInputType *MaterialModelInputs<dim>::get_additional_input() const
1477  {
1478  for (unsigned int i=0; i<additional_inputs.size(); ++i)
1479  {
1480  const AdditionalInputType *result = dynamic_cast<const AdditionalInputType *> (additional_inputs[i].get());
1481  if (result)
1482  return result;
1483  }
1484  return nullptr;
1485  }
1486 
1487 
1488  template <int dim>
1489  template <class AdditionalOutputType>
1491  {
1492  for (unsigned int i=0; i<additional_outputs.size(); ++i)
1493  {
1494  AdditionalOutputType *result = dynamic_cast<AdditionalOutputType *> (additional_outputs[i].get());
1495  if (result)
1496  return result;
1497  }
1498  return nullptr;
1499  }
1500 
1501 
1502  template <int dim>
1503  template <class AdditionalOutputType>
1504  const AdditionalOutputType *MaterialModelOutputs<dim>::get_additional_output() const
1505  {
1506  for (unsigned int i=0; i<additional_outputs.size(); ++i)
1507  {
1508  const AdditionalOutputType *result = dynamic_cast<const AdditionalOutputType *> (additional_outputs[i].get());
1509  if (result)
1510  return result;
1511  }
1512  return nullptr;
1513  }
1514 
1515 
1516  template <int dim>
1518  {
1519  Assert(this->additional_outputs.empty(), ExcMessage("Destination of move needs to be empty!"));
1520  this->additional_outputs = std::move(other.additional_outputs);
1521  }
1522 
1523 
1531 #define ASPECT_REGISTER_MATERIAL_MODEL(classname,name,description) \
1532  template class classname<2>; \
1533  template class classname<3>; \
1534  namespace ASPECT_REGISTER_MATERIAL_MODEL_ ## classname \
1535  { \
1536  aspect::internal::Plugins::RegisterHelper<aspect::MaterialModel::Interface<2>,classname<2>> \
1537  dummy_ ## classname ## _2d (&aspect::MaterialModel::register_material_model<2>, \
1538  name, description); \
1539  aspect::internal::Plugins::RegisterHelper<aspect::MaterialModel::Interface<3>,classname<3>> \
1540  dummy_ ## classname ## _3d (&aspect::MaterialModel::register_material_model<3>, \
1541  name, description); \
1542  }
1543  }
1544 }
1545 
1546 
1547 #endif
void write_plugin_graph(std::ostream &output_stream)
void average_property(const AveragingOperation operation, const FullMatrix< double > &projection_matrix, const FullMatrix< double > &expansion_matrix, std::vector< double > &values_out)
::TrilinosWrappers::MPI::BlockVector BlockVector
Definition: global.h:229
AdditionalInputType * get_additional_input()
Definition: interface.h:1462
std::vector< SymmetricTensor< 2, dim > > strain_rate
Definition: interface.h:411
std::vector< double > thermal_expansion_coefficients
Definition: interface.h:539
virtual void average(const MaterialAveraging::AveragingOperation, const FullMatrix< double > &, const FullMatrix< double > &)
Definition: interface.h:820
MaterialProperties::Property requested_properties
Definition: interface.h:437
std::vector< std::vector< double > > prescribed_field_outputs
Definition: interface.h:1027
std::vector< std::vector< double > > composition
Definition: interface.h:397
std::vector< Point< dim > > position
Definition: interface.h:365
ElasticOutputs(const unsigned int n_points)
Definition: interface.h:1155
AdditionalOutputType * get_additional_output()
Definition: interface.h:1490
std::vector< std::vector< double > > reaction_rates
Definition: interface.h:974
void declare_parameters(ParameterHandler &prm)
bool identifies_single_variable(const Dependence dependence)
void average(const MaterialAveraging::AveragingOperation, const FullMatrix< double > &, const FullMatrix< double > &) override
Definition: interface.h:1078
void average(const AveragingOperation operation, const typename DoFHandler< dim >::active_cell_iterator &cell, const Quadrature< dim > &quadrature_formula, const Mapping< dim > &mapping, MaterialModelOutputs< dim > &values_out)
const std::vector< std::string > names
Definition: interface.h:905
Property operator|=(Property &d1, const Property d2)
Definition: interface.h:229
std::vector< std::vector< double > > reaction_terms
Definition: interface.h:608
AveragingOperation parse_averaging_operation_name(const std::string &s)
DoFHandler< dim >::active_cell_iterator current_cell
Definition: interface.h:428
Property operator|(const Property d1, const Property d2)
Definition: interface.h:223
AdditionalMaterialOutputsStokesRHS(const unsigned int n_points)
Definition: interface.h:1071
std::vector< std::unique_ptr< AdditionalMaterialInputs< dim > > > additional_inputs
Definition: interface.h:444
void move_additional_outputs_from(MaterialModelOutputs< dim > &other)
Definition: interface.h:1517
std::vector< Tensor< 1, dim > > pressure_gradient
Definition: interface.h:381
std::vector< std::vector< double > > output_values
Definition: interface.h:902
Interface< dim > * create_material_model(const std::string &model_name)
std::vector< Tensor< 2, dim > > elastic_force
Definition: interface.h:1175
void average(const MaterialAveraging::AveragingOperation operation, const FullMatrix< double > &, const FullMatrix< double > &) override
Definition: interface.h:1162
std::vector< double > entropy_derivative_pressure
Definition: interface.h:562
std::vector< double > thermal_conductivities
Definition: interface.h:549
std::vector< double > entropy_derivative_temperature
Definition: interface.h:570
NonlinearDependence::ModelDependence model_dependence
Definition: interface.h:1365
std::string get_valid_model_names_pattern()
std::vector< Tensor< 1, dim > > velocity
Definition: interface.h:390
void register_material_model(const std::string &name, const std::string &description, void(*declare_parameters_function)(ParameterHandler &), Interface< dim > *(*factory_function)())
void average(const MaterialAveraging::AveragingOperation, const FullMatrix< double > &, const FullMatrix< double > &) override
Definition: interface.h:892
std::vector< std::unique_ptr< AdditionalMaterialOutputs< dim > > > additional_outputs
Definition: interface.h:615