ASPECT
interface.h
Go to the documentation of this file.
1 /*
2  Copyright (C) 2011 - 2018 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 #include <deal.II/base/quadrature.h>
31 #include <deal.II/base/parameter_handler.h>
32 #include <deal.II/dofs/dof_handler.h>
33 #include <deal.II/dofs/dof_accessor.h>
34 #include <deal.II/fe/mapping.h>
35 #include <deal.II/fe/fe_values.h>
36 #include <deal.II/fe/component_mask.h>
37 #include <deal.II/numerics/data_postprocessor.h>
38 #include <deal.II/base/signaling_nan.h>
39 
40 namespace aspect
41 {
42  template <int dim>
43  struct Introspection;
44 
45 
53  namespace MaterialModel
54  {
55  using namespace dealii;
56 
61  namespace NonlinearDependence
62  {
85  {
87 
88  none = 1,
90  pressure = 4,
93 
95  };
96 
97 
102  const Dependence d2)
103  {
104  return Dependence(static_cast<int>(d1) | static_cast<int>(d2));
105  }
106 
108  const Dependence d2)
109  {
110  d1 = (d1 | d2);
111  return d1;
112  }
113 
119  {
125 
131 
137 
143 
149 
155  ModelDependence ();
156  };
157 
166  bool
167  identifies_single_variable(const Dependence dependence);
168  }
169 
170 
171  template <int dim> class AdditionalMaterialInputs;
172 
173 
180  template <int dim>
182  {
193  MaterialModelInputs(const unsigned int n_points,
194  const unsigned int n_comp);
195 
205  const Introspection<dim> &introspection,
206  const bool use_strain_rate = true);
207 
208 
222  const typename DoFHandler<dim>::active_cell_iterator &cell,
223  const Introspection<dim> &introspection,
224  const LinearAlgebra::BlockVector &solution_vector,
225  const bool use_strain_rates = true);
226 
241 
246 
251  MaterialModelInputs &operator= (const MaterialModelInputs &source) = delete;
252 
256  MaterialModelInputs &operator= (MaterialModelInputs &&) = default;
257 
262  void reinit(const FEValuesBase<dim,dim> &fe_values,
263  const typename DoFHandler<dim>::active_cell_iterator &cell,
264  const Introspection<dim> &introspection,
265  const LinearAlgebra::BlockVector &solution_vector,
266  const bool use_strain_rates = true);
267 
268 
273  std::vector<Point<dim> > position;
274 
278  std::vector<double> temperature;
279 
283  std::vector<double> pressure;
284 
289  std::vector<Tensor<1,dim> > pressure_gradient;
290 
298  std::vector<Tensor<1,dim> > velocity;
299 
305  std::vector<std::vector<double> > composition;
306 
319  std::vector<SymmetricTensor<2,dim> > strain_rate;
320 
334 
351 
357  std::vector<std::unique_ptr<AdditionalMaterialInputs<dim> > > additional_inputs;
358 
365  template <class AdditionalInputType>
366  AdditionalInputType *get_additional_input();
367 
371  template <class AdditionalInputType>
372  const AdditionalInputType *get_additional_input() const;
373  };
374 
375 
376  template <int dim> class AdditionalMaterialOutputs;
377 
378 
385  template <int dim>
387  {
398  MaterialModelOutputs (const unsigned int n_points,
399  const unsigned int n_comp);
400 
415 
420 
425  MaterialModelOutputs &operator= (const MaterialModelOutputs &source) = delete;
426 
430  MaterialModelOutputs &operator= (MaterialModelOutputs &&) = default;
431 
432 
436  std::vector<double> viscosities;
437 
441  std::vector<double> densities;
442 
447  std::vector<double> thermal_expansion_coefficients;
448 
452  std::vector<double> specific_heat;
453 
457  std::vector<double> thermal_conductivities;
458 
463  std::vector<double> compressibilities;
464 
470  std::vector<double> entropy_derivative_pressure;
471 
478  std::vector<double> entropy_derivative_temperature;
479 
516  std::vector<std::vector<double> > reaction_terms;
517 
523  std::vector<std::unique_ptr<AdditionalMaterialOutputs<dim> > > additional_outputs;
524 
532  template <class AdditionalOutputType>
533  AdditionalOutputType *get_additional_output();
534 
538  template <class AdditionalOutputType>
539  const AdditionalOutputType *get_additional_output() const;
540 
545  void move_additional_outputs_from(MaterialModelOutputs<dim> &other);
546  };
547 
548 
549 
564  namespace MaterialAveraging
565  {
604  {
613  };
614 
615 
622  std::string get_averaging_operation_names ();
623 
629  AveragingOperation parse_averaging_operation_name (const std::string &s);
630 
636  template <int dim>
637  void average (const AveragingOperation operation,
638  const typename DoFHandler<dim>::active_cell_iterator &cell,
639  const Quadrature<dim> &quadrature_formula,
640  const Mapping<dim> &mapping,
641  MaterialModelOutputs<dim> &values_out);
642 
648  void average_property (const AveragingOperation operation,
649  const FullMatrix<double> &projection_matrix,
650  const FullMatrix<double> &expansion_matrix,
651  std::vector<double> &values_out);
652  }
653 
654 
673  template <int dim>
675  {
676  public:
678  {}
679 
684  virtual void
685  fill (const LinearAlgebra::BlockVector &solution,
686  const FEValuesBase<dim> &fe_values,
687  const Introspection<dim> &introspection) = 0;
688  };
689 
690 
717  template <int dim>
719  {
720  public:
725  virtual ~AdditionalMaterialOutputs() = default;
726 
727  virtual void average (const MaterialAveraging::AveragingOperation /*operation*/,
728  const FullMatrix<double> &/*projection_matrix*/,
729  const FullMatrix<double> &/*expansion_matrix*/)
730  {}
731  };
732 
733 
755  template <int dim>
757  {
758  public:
767  NamedAdditionalMaterialOutputs(const std::vector<std::string> &output_names);
768 
772  ~NamedAdditionalMaterialOutputs() override;
773 
778  const std::vector<std::string> &get_names() const;
779 
784  virtual std::vector<double> get_nth_output(const unsigned int idx) const = 0;
785 
787  const FullMatrix<double> &/*projection_matrix*/,
788  const FullMatrix<double> &/*expansion_matrix*/) override
789  {}
790 
791  private:
792  const std::vector<std::string> names;
793  };
794 
795 
801  template <int dim>
803  {
804  public:
805  SeismicAdditionalOutputs(const unsigned int n_points);
806 
807  std::vector<double> get_nth_output(const unsigned int idx) const override;
808 
814  std::vector<double> vs;
815 
821  std::vector<double> vp;
822  };
823 
824 
845  template <int dim>
847  {
848  public:
849  ReactionRateOutputs (const unsigned int n_points,
850  const unsigned int n_comp);
851 
852  std::vector<double> get_nth_output(const unsigned int idx) const override;
853 
861  std::vector<std::vector<double> > reaction_rates;
862  };
863 
864 
883  template <int dim>
885  {
886  public:
887  PrescribedFieldOutputs (const unsigned int n_points,
888  const unsigned int n_comp);
889 
890  std::vector<double> get_nth_output(const unsigned int idx) const override;
891 
899  std::vector<std::vector<double> > prescribed_field_outputs;
900  };
901 
915  template <int dim>
917  {
918  public:
919  PrescribedTemperatureOutputs (const unsigned int n_points);
920 
921  std::vector<double> get_nth_output(const unsigned int idx) const override;
922 
930  std::vector<double> prescribed_temperature_outputs;
931  };
932 
939  template <int dim>
941  {
942  public:
943  AdditionalMaterialOutputsStokesRHS(const unsigned int n_points)
944  : rhs_u(n_points), rhs_p(n_points), rhs_melt_pc(n_points)
945  {}
946 
948  {}
949 
951  const FullMatrix<double> &/*projection_matrix*/,
952  const FullMatrix<double> &/*expansion_matrix*/) override
953  {
954  // TODO: not implemented
955  }
956 
962  std::vector<Tensor<1,dim> > rhs_u;
963 
968  std::vector<double> rhs_p;
969 
974  std::vector<double> rhs_melt_pc;
975  };
976 
977 
978 
994  template <int dim>
996  {
997  public:
1001  explicit PrescribedPlasticDilation (const unsigned int n_points);
1002 
1006  virtual std::vector<double> get_nth_output(const unsigned int idx) const;
1007 
1012  std::vector<double> dilation;
1013  };
1014 
1015 
1016 
1023  template <int dim>
1025  {
1026  public:
1027  ElasticOutputs(const unsigned int n_points)
1028  : elastic_force(n_points, numbers::signaling_nan<Tensor<2,dim> >() )
1029  {}
1030 
1031  ~ElasticOutputs() override
1032  {}
1033 
1035  const FullMatrix<double> &/*projection_matrix*/,
1036  const FullMatrix<double> &/*expansion_matrix*/) override
1037  {
1039  return;
1040  }
1041 
1047  std::vector<Tensor<2,dim> > elastic_force;
1048  };
1049 
1050 
1051 
1067  template <int dim>
1068  class Interface
1069  {
1070  public:
1085 
1090  virtual ~Interface();
1091 
1097  virtual
1098  void
1099  initialize ();
1100 
1105  virtual void update ();
1106 
1119  get_model_dependence () const;
1120 
1130  virtual bool is_compressible () const = 0;
1157  virtual double reference_viscosity () const = 0;
1167  virtual
1168  void evaluate (const MaterialModel::MaterialModelInputs<dim> &in,
1180  static
1181  void
1183 
1190  virtual
1191  void
1192  parse_parameters (ParameterHandler &prm);
1202  virtual
1203  void
1204  create_additional_named_outputs (MaterialModelOutputs &outputs) const;
1205 
1206 
1214  virtual
1215  void
1216  fill_additional_material_model_inputs(MaterialModel::MaterialModelInputs<dim> &input,
1217  const LinearAlgebra::BlockVector &solution,
1218  const FEValuesBase<dim> &fe_values,
1219  const Introspection<dim> &introspection) const;
1220 
1221  protected:
1238  };
1239 
1255  template <int dim>
1256  void
1257  register_material_model (const std::string &name,
1258  const std::string &description,
1259  void (*declare_parameters_function) (ParameterHandler &),
1260  Interface<dim> *(*factory_function) ());
1261 
1272  template <int dim>
1273  Interface<dim> *
1274  create_material_model (const std::string &model_name);
1275 
1276 
1288  template <int dim>
1289  Interface<dim> *
1291 
1292 
1299  template <int dim>
1300  std::string
1302 
1303 
1309  template <int dim>
1310  void
1312 
1313 
1314 
1324  template <int dim>
1325  void
1326  write_plugin_graph (std::ostream &output_stream);
1327 
1328 
1329 
1330 // --------------------- template function definitions ----------------------------------
1331 
1332  template <int dim>
1333  template <class AdditionalInputType>
1335  {
1336  for (unsigned int i=0; i<additional_inputs.size(); ++i)
1337  {
1338  AdditionalInputType *result = dynamic_cast<AdditionalInputType *> (additional_inputs[i].get());
1339  if (result)
1340  return result;
1341  }
1342  return nullptr;
1343  }
1344 
1345 
1346  template <int dim>
1347  template <class AdditionalInputType>
1348  const AdditionalInputType *MaterialModelInputs<dim>::get_additional_input() const
1349  {
1350  for (unsigned int i=0; i<additional_inputs.size(); ++i)
1351  {
1352  const AdditionalInputType *result = dynamic_cast<const AdditionalInputType *> (additional_inputs[i].get());
1353  if (result)
1354  return result;
1355  }
1356  return nullptr;
1357  }
1358 
1359 
1360  template <int dim>
1361  template <class AdditionalOutputType>
1363  {
1364  for (unsigned int i=0; i<additional_outputs.size(); ++i)
1365  {
1366  AdditionalOutputType *result = dynamic_cast<AdditionalOutputType *> (additional_outputs[i].get());
1367  if (result)
1368  return result;
1369  }
1370  return nullptr;
1371  }
1372 
1373 
1374  template <int dim>
1375  template <class AdditionalOutputType>
1376  const AdditionalOutputType *MaterialModelOutputs<dim>::get_additional_output() const
1377  {
1378  for (unsigned int i=0; i<additional_outputs.size(); ++i)
1379  {
1380  const AdditionalOutputType *result = dynamic_cast<const AdditionalOutputType *> (additional_outputs[i].get());
1381  if (result)
1382  return result;
1383  }
1384  return nullptr;
1385  }
1386 
1387 
1388  template <int dim>
1390  {
1391  Assert(this->additional_outputs.empty(), ExcMessage("Destination of move needs to be empty!"));
1392  this->additional_outputs = std::move(other.additional_outputs);
1393  }
1394 
1395 
1403 #define ASPECT_REGISTER_MATERIAL_MODEL(classname,name,description) \
1404  template class classname<2>; \
1405  template class classname<3>; \
1406  namespace ASPECT_REGISTER_MATERIAL_MODEL_ ## classname \
1407  { \
1408  aspect::internal::Plugins::RegisterHelper<aspect::MaterialModel::Interface<2>,classname<2> > \
1409  dummy_ ## classname ## _2d (&aspect::MaterialModel::register_material_model<2>, \
1410  name, description); \
1411  aspect::internal::Plugins::RegisterHelper<aspect::MaterialModel::Interface<3>,classname<3> > \
1412  dummy_ ## classname ## _3d (&aspect::MaterialModel::register_material_model<3>, \
1413  name, description); \
1414  }
1415  }
1416 }
1417 
1418 
1419 #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)
AdditionalInputType * get_additional_input()
Definition: interface.h:1334
std::vector< SymmetricTensor< 2, dim > > strain_rate
Definition: interface.h:319
std::vector< double > thermal_expansion_coefficients
Definition: interface.h:447
virtual void average(const MaterialAveraging::AveragingOperation, const FullMatrix< double > &, const FullMatrix< double > &)
Definition: interface.h:727
std::vector< std::vector< double > > prescribed_field_outputs
Definition: interface.h:899
MaterialModel::MaterialModelOutputs< dim > MaterialModelOutputs
Definition: interface.h:1084
std::vector< std::vector< double > > composition
Definition: interface.h:305
std::vector< Point< dim > > position
Definition: interface.h:273
#define AssertThrow(cond, exc)
ElasticOutputs(const unsigned int n_points)
Definition: interface.h:1027
AdditionalOutputType * get_additional_output()
Definition: interface.h:1362
std::vector< std::vector< double > > reaction_rates
Definition: interface.h:861
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:950
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)
UpdateFlags operator|(const UpdateFlags f1, const UpdateFlags f2)
const std::vector< std::string > names
Definition: interface.h:792
static ::ExceptionBase & ExcMessage(std::string arg1)
std::vector< std::vector< double > > reaction_terms
Definition: interface.h:516
MaterialModel::MaterialModelInputs< dim > MaterialModelInputs
Definition: interface.h:1077
AveragingOperation parse_averaging_operation_name(const std::string &s)
DoFHandler< dim >::active_cell_iterator current_cell
Definition: interface.h:350
AdditionalMaterialOutputsStokesRHS(const unsigned int n_points)
Definition: interface.h:943
#define Assert(cond, exc)
std::vector< std::unique_ptr< AdditionalMaterialInputs< dim > > > additional_inputs
Definition: interface.h:357
void move_additional_outputs_from(MaterialModelOutputs< dim > &other)
Definition: interface.h:1389
std::vector< Tensor< 1, dim > > pressure_gradient
Definition: interface.h:289
Interface< dim > * create_material_model(const std::string &model_name)
std::vector< Tensor< 2, dim > > elastic_force
Definition: interface.h:1047
void average(const MaterialAveraging::AveragingOperation operation, const FullMatrix< double > &, const FullMatrix< double > &) override
Definition: interface.h:1034
Definition: compat.h:53
std::vector< double > entropy_derivative_pressure
Definition: interface.h:470
UpdateFlags & operator|=(UpdateFlags &f1, const UpdateFlags f2)
std::vector< double > thermal_conductivities
Definition: interface.h:457
std::vector< double > entropy_derivative_temperature
Definition: interface.h:478
NonlinearDependence::ModelDependence model_dependence
Definition: interface.h:1237
std::string get_valid_model_names_pattern()
static ::ExceptionBase & ExcNotImplemented()
std::vector< Tensor< 1, dim > > velocity
Definition: interface.h:298
void register_material_model(const std::string &name, const std::string &description, void(*declare_parameters_function)(ParameterHandler &), Interface< dim > *(*factory_function)())
const DoFHandler< dim >::active_cell_iterator *cell DEAL_II_DEPRECATED
Definition: interface.h:333
void average(const MaterialAveraging::AveragingOperation, const FullMatrix< double > &, const FullMatrix< double > &) override
Definition: interface.h:786
std::vector< std::unique_ptr< AdditionalMaterialOutputs< dim > > > additional_outputs
Definition: interface.h:523