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 #if !DEAL_II_VERSION_GTE(9,3,0)
37 #endif
39 
44 #include <deal.II/fe/mapping.h>
45 #include <deal.II/fe/fe_values.h>
49 
50 namespace aspect
51 {
52  template <int dim>
53  struct Introspection;
54 
55 
63  namespace MaterialModel
64  {
65  using namespace dealii;
66 
71  namespace NonlinearDependence
72  {
95  {
97 
98  none = 1,
100  pressure = 4,
103 
105  };
106 
107 
112  const Dependence d2)
113  {
114  return Dependence(static_cast<int>(d1) | static_cast<int>(d2));
115  }
116 
118  const Dependence d2)
119  {
120  d1 = (d1 | d2);
121  return d1;
122  }
123 
129  {
135 
141 
147 
153 
159 
165  ModelDependence ();
166  };
167 
176  bool
177  identifies_single_variable(const Dependence dependence);
178  }
179 
184  namespace MaterialProperties
185  {
195  enum Property
196  {
198 
199  none = 1,
201  density = 4,
209 
212  specific_heat |
217  viscosity |
219  };
220 
225  inline Property operator | (const Property d1,
226  const Property d2)
227  {
228  return Property(static_cast<int>(d1) | static_cast<int>(d2));
229  }
230 
232  const Property d2)
233  {
234  return (d1 | d2);
235  }
236  }
237 
238  template <int dim> class AdditionalMaterialInputs;
239 
246  template <int dim>
248  {
259  MaterialModelInputs(const unsigned int n_points,
260  const unsigned int n_comp);
261 
271  const Introspection<dim> &introspection,
272  const bool use_strain_rate = true);
273 
274 
288  const typename DoFHandler<dim>::active_cell_iterator &cell,
289  const Introspection<dim> &introspection,
290  const LinearAlgebra::BlockVector &solution_vector,
291  const bool use_strain_rates = true);
292 
307 
312 
317  MaterialModelInputs &operator= (const MaterialModelInputs &source) = delete;
318 
322  MaterialModelInputs &operator= (MaterialModelInputs &&) = default;
323 
328  void reinit(const FEValuesBase<dim,dim> &fe_values,
329  const typename DoFHandler<dim>::active_cell_iterator &cell,
330  const Introspection<dim> &introspection,
331  const LinearAlgebra::BlockVector &solution_vector,
332  const bool use_strain_rates = true);
333 
338  unsigned int n_evaluation_points() const;
339 
346  bool requests_property(const MaterialProperties::Property &property) const;
347 
352  std::vector<Point<dim> > position;
353 
357  std::vector<double> temperature;
358 
362  std::vector<double> pressure;
363 
368  std::vector<Tensor<1,dim> > pressure_gradient;
369 
377  std::vector<Tensor<1,dim> > velocity;
378 
384  std::vector<std::vector<double> > composition;
385 
398  std::vector<SymmetricTensor<2,dim> > strain_rate;
399 
416 
425 
431  std::vector<std::unique_ptr<AdditionalMaterialInputs<dim> > > additional_inputs;
432 
439  template <class AdditionalInputType>
440  AdditionalInputType *get_additional_input();
441 
445  template <class AdditionalInputType>
446  const AdditionalInputType *get_additional_input() const;
447  };
448 
449 
450  template <int dim> class AdditionalMaterialOutputs;
451 
452 
459  template <int dim>
461  {
472  MaterialModelOutputs (const unsigned int n_points,
473  const unsigned int n_comp);
474 
489 
494 
499  MaterialModelOutputs &operator= (const MaterialModelOutputs &source) = delete;
500 
504  MaterialModelOutputs &operator= (MaterialModelOutputs &&) = default;
505 
510  unsigned int n_evaluation_points() const;
511 
515  std::vector<double> viscosities;
516 
520  std::vector<double> densities;
521 
526  std::vector<double> thermal_expansion_coefficients;
527 
531  std::vector<double> specific_heat;
532 
536  std::vector<double> thermal_conductivities;
537 
542  std::vector<double> compressibilities;
543 
549  std::vector<double> entropy_derivative_pressure;
550 
557  std::vector<double> entropy_derivative_temperature;
558 
595  std::vector<std::vector<double> > reaction_terms;
596 
602  std::vector<std::unique_ptr<AdditionalMaterialOutputs<dim> > > additional_outputs;
603 
611  template <class AdditionalOutputType>
612  AdditionalOutputType *get_additional_output();
613 
617  template <class AdditionalOutputType>
618  const AdditionalOutputType *get_additional_output() const;
619 
624  void move_additional_outputs_from(MaterialModelOutputs<dim> &other);
625  };
626 
627 
628 
643  namespace MaterialAveraging
644  {
684  {
694  };
695 
696 
703  std::string get_averaging_operation_names ();
704 
710  AveragingOperation parse_averaging_operation_name (const std::string &s);
711 
717  template <int dim>
718  void average (const AveragingOperation operation,
719  const typename DoFHandler<dim>::active_cell_iterator &cell,
720  const Quadrature<dim> &quadrature_formula,
721  const Mapping<dim> &mapping,
722  MaterialModelOutputs<dim> &values_out);
723 
729  void average_property (const AveragingOperation operation,
730  const FullMatrix<double> &projection_matrix,
731  const FullMatrix<double> &expansion_matrix,
732  std::vector<double> &values_out);
733  }
734 
735 
754  template <int dim>
756  {
757  public:
759  {}
760 
765  virtual void
766  fill (const LinearAlgebra::BlockVector &solution,
767  const FEValuesBase<dim> &fe_values,
768  const Introspection<dim> &introspection) = 0;
769  };
770 
771 
798  template <int dim>
800  {
801  public:
806  virtual ~AdditionalMaterialOutputs() = default;
807 
808  virtual void average (const MaterialAveraging::AveragingOperation /*operation*/,
809  const FullMatrix<double> &/*projection_matrix*/,
810  const FullMatrix<double> &/*expansion_matrix*/)
811  {}
812  };
813 
814 
836  template <int dim>
838  {
839  public:
848  NamedAdditionalMaterialOutputs(const std::vector<std::string> &output_names);
849 
860  NamedAdditionalMaterialOutputs(const std::vector<std::string> &output_names,
861  const unsigned int n_points);
862 
866  ~NamedAdditionalMaterialOutputs() override;
867 
872  const std::vector<std::string> &get_names() const;
873 
878  virtual std::vector<double> get_nth_output(const unsigned int idx) const;
879 
881  const FullMatrix<double> &/*projection_matrix*/,
882  const FullMatrix<double> &/*expansion_matrix*/) override
883  {}
884 
885 
890  std::vector<std::vector<double> > output_values;
891 
892  private:
893  const std::vector<std::string> names;
894  };
895 
896 
902  template <int dim>
904  {
905  public:
906  SeismicAdditionalOutputs(const unsigned int n_points);
907 
908  std::vector<double> get_nth_output(const unsigned int idx) const override;
909 
915  std::vector<double> vs;
916 
922  std::vector<double> vp;
923  };
924 
925 
946  template <int dim>
948  {
949  public:
950  ReactionRateOutputs (const unsigned int n_points,
951  const unsigned int n_comp);
952 
953  std::vector<double> get_nth_output(const unsigned int idx) const override;
954 
962  std::vector<std::vector<double> > reaction_rates;
963  };
964 
965 
984  template <int dim>
986  {
987  public:
988  PrescribedFieldOutputs (const unsigned int n_points,
989  const unsigned int n_comp);
990 
991  std::vector<double> get_nth_output(const unsigned int idx) const override;
992 
1000  std::vector<std::vector<double> > prescribed_field_outputs;
1001  };
1002 
1016  template <int dim>
1018  {
1019  public:
1020  PrescribedTemperatureOutputs (const unsigned int n_points);
1021 
1022  std::vector<double> get_nth_output(const unsigned int idx) const override;
1023 
1031  std::vector<double> prescribed_temperature_outputs;
1032  };
1033 
1040  template <int dim>
1042  {
1043  public:
1044  AdditionalMaterialOutputsStokesRHS(const unsigned int n_points)
1045  : rhs_u(n_points), rhs_p(n_points), rhs_melt_pc(n_points)
1046  {}
1047 
1049  {}
1050 
1052  const FullMatrix<double> &/*projection_matrix*/,
1053  const FullMatrix<double> &/*expansion_matrix*/) override
1054  {
1055  // TODO: not implemented
1056  }
1057 
1063  std::vector<Tensor<1,dim> > rhs_u;
1064 
1069  std::vector<double> rhs_p;
1070 
1075  std::vector<double> rhs_melt_pc;
1076  };
1077 
1078 
1079 
1095  template <int dim>
1097  {
1098  public:
1102  explicit PrescribedPlasticDilation (const unsigned int n_points);
1103 
1107  virtual std::vector<double> get_nth_output(const unsigned int idx) const;
1108 
1113  std::vector<double> dilation;
1114  };
1115 
1116 
1117 
1124  template <int dim>
1126  {
1127  public:
1128  ElasticOutputs(const unsigned int n_points)
1129  : elastic_force(n_points, numbers::signaling_nan<Tensor<2,dim> >() )
1130  {}
1131 
1132  ~ElasticOutputs() override
1133  {}
1134 
1136  const FullMatrix<double> &/*projection_matrix*/,
1137  const FullMatrix<double> &/*expansion_matrix*/) override
1138  {
1140  return;
1141  }
1142 
1148  std::vector<Tensor<2,dim> > elastic_force;
1149  };
1150 
1151 
1152 
1168  template <int dim>
1169  class Interface
1170  {
1171  public:
1186 
1191  virtual ~Interface();
1192 
1198  virtual
1199  void
1200  initialize ();
1201 
1206  virtual void update ();
1207 
1220  get_model_dependence () const;
1221 
1231  virtual bool is_compressible () const = 0;
1258  virtual double reference_viscosity () const = 0;
1268  virtual
1269  void evaluate (const MaterialModel::MaterialModelInputs<dim> &in,
1281  static
1282  void
1284 
1291  virtual
1292  void
1293  parse_parameters (ParameterHandler &prm);
1303  virtual
1304  void
1305  create_additional_named_outputs (MaterialModelOutputs &outputs) const;
1306 
1307 
1315  virtual
1316  void
1317  fill_additional_material_model_inputs(MaterialModel::MaterialModelInputs<dim> &input,
1318  const LinearAlgebra::BlockVector &solution,
1319  const FEValuesBase<dim> &fe_values,
1320  const Introspection<dim> &introspection) const;
1321 
1322  protected:
1339  };
1340 
1356  template <int dim>
1357  void
1358  register_material_model (const std::string &name,
1359  const std::string &description,
1360  void (*declare_parameters_function) (ParameterHandler &),
1361  Interface<dim> *(*factory_function) ());
1362 
1373  template <int dim>
1374  Interface<dim> *
1375  create_material_model (const std::string &model_name);
1376 
1377 
1389  template <int dim>
1390  Interface<dim> *
1392 
1393 
1400  template <int dim>
1401  std::string
1403 
1404 
1410  template <int dim>
1411  void
1413 
1414 
1415 
1425  template <int dim>
1426  void
1427  write_plugin_graph (std::ostream &output_stream);
1428 
1429 
1430 
1431 // --------------------- template function definitions ----------------------------------
1432 
1433  template <int dim>
1434  template <class AdditionalInputType>
1436  {
1437  for (unsigned int i=0; i<additional_inputs.size(); ++i)
1438  {
1439  AdditionalInputType *result = dynamic_cast<AdditionalInputType *> (additional_inputs[i].get());
1440  if (result)
1441  return result;
1442  }
1443  return nullptr;
1444  }
1445 
1446 
1447  template <int dim>
1448  template <class AdditionalInputType>
1449  const AdditionalInputType *MaterialModelInputs<dim>::get_additional_input() const
1450  {
1451  for (unsigned int i=0; i<additional_inputs.size(); ++i)
1452  {
1453  const AdditionalInputType *result = dynamic_cast<const AdditionalInputType *> (additional_inputs[i].get());
1454  if (result)
1455  return result;
1456  }
1457  return nullptr;
1458  }
1459 
1460 
1461  template <int dim>
1462  template <class AdditionalOutputType>
1464  {
1465  for (unsigned int i=0; i<additional_outputs.size(); ++i)
1466  {
1467  AdditionalOutputType *result = dynamic_cast<AdditionalOutputType *> (additional_outputs[i].get());
1468  if (result)
1469  return result;
1470  }
1471  return nullptr;
1472  }
1473 
1474 
1475  template <int dim>
1476  template <class AdditionalOutputType>
1477  const AdditionalOutputType *MaterialModelOutputs<dim>::get_additional_output() const
1478  {
1479  for (unsigned int i=0; i<additional_outputs.size(); ++i)
1480  {
1481  const AdditionalOutputType *result = dynamic_cast<const AdditionalOutputType *> (additional_outputs[i].get());
1482  if (result)
1483  return result;
1484  }
1485  return nullptr;
1486  }
1487 
1488 
1489  template <int dim>
1491  {
1492  Assert(this->additional_outputs.empty(), ExcMessage("Destination of move needs to be empty!"));
1493  this->additional_outputs = std::move(other.additional_outputs);
1494  }
1495 
1496 
1504 #define ASPECT_REGISTER_MATERIAL_MODEL(classname,name,description) \
1505  template class classname<2>; \
1506  template class classname<3>; \
1507  namespace ASPECT_REGISTER_MATERIAL_MODEL_ ## classname \
1508  { \
1509  aspect::internal::Plugins::RegisterHelper<aspect::MaterialModel::Interface<2>,classname<2> > \
1510  dummy_ ## classname ## _2d (&aspect::MaterialModel::register_material_model<2>, \
1511  name, description); \
1512  aspect::internal::Plugins::RegisterHelper<aspect::MaterialModel::Interface<3>,classname<3> > \
1513  dummy_ ## classname ## _3d (&aspect::MaterialModel::register_material_model<3>, \
1514  name, description); \
1515  }
1516  }
1517 }
1518 
1519 
1520 #endif
void write_plugin_graph(std::ostream &output_stream)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
void average_property(const AveragingOperation operation, const FullMatrix< double > &projection_matrix, const FullMatrix< double > &expansion_matrix, std::vector< double > &values_out)
AdditionalInputType * get_additional_input()
Definition: interface.h:1435
std::vector< SymmetricTensor< 2, dim > > strain_rate
Definition: interface.h:398
std::vector< double > thermal_expansion_coefficients
Definition: interface.h:526
virtual void average(const MaterialAveraging::AveragingOperation, const FullMatrix< double > &, const FullMatrix< double > &)
Definition: interface.h:808
MaterialProperties::Property requested_properties
Definition: interface.h:424
std::vector< std::vector< double > > prescribed_field_outputs
Definition: interface.h:1000
std::vector< std::vector< double > > composition
Definition: interface.h:384
std::vector< Point< dim > > position
Definition: interface.h:352
#define AssertThrow(cond, exc)
ElasticOutputs(const unsigned int n_points)
Definition: interface.h:1128
AdditionalOutputType * get_additional_output()
Definition: interface.h:1463
std::vector< std::vector< double > > reaction_rates
Definition: interface.h:962
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:1051
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:893
Property operator|=(Property &d1, const Property d2)
Definition: interface.h:231
static ::ExceptionBase & ExcMessage(std::string arg1)
std::vector< std::vector< double > > reaction_terms
Definition: interface.h:595
AveragingOperation parse_averaging_operation_name(const std::string &s)
DoFHandler< dim >::active_cell_iterator current_cell
Definition: interface.h:415
Property operator|(const Property d1, const Property d2)
Definition: interface.h:225
AdditionalMaterialOutputsStokesRHS(const unsigned int n_points)
Definition: interface.h:1044
#define Assert(cond, exc)
std::vector< std::unique_ptr< AdditionalMaterialInputs< dim > > > additional_inputs
Definition: interface.h:431
void move_additional_outputs_from(MaterialModelOutputs< dim > &other)
Definition: interface.h:1490
std::vector< Tensor< 1, dim > > pressure_gradient
Definition: interface.h:368
std::vector< std::vector< double > > output_values
Definition: interface.h:890
Interface< dim > * create_material_model(const std::string &model_name)
std::vector< Tensor< 2, dim > > elastic_force
Definition: interface.h:1148
void average(const MaterialAveraging::AveragingOperation operation, const FullMatrix< double > &, const FullMatrix< double > &) override
Definition: interface.h:1135
std::vector< double > entropy_derivative_pressure
Definition: interface.h:549
std::vector< double > thermal_conductivities
Definition: interface.h:536
std::vector< double > entropy_derivative_temperature
Definition: interface.h:557
T signaling_nan()
NonlinearDependence::ModelDependence model_dependence
Definition: interface.h:1338
std::string get_valid_model_names_pattern()
static ::ExceptionBase & ExcNotImplemented()
std::vector< Tensor< 1, dim > > velocity
Definition: interface.h:377
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:880
std::vector< std::unique_ptr< AdditionalMaterialOutputs< dim > > > additional_outputs
Definition: interface.h:602