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  };
541 
542 
543 
558  namespace MaterialAveraging
559  {
595  {
603  };
604 
605 
612  std::string get_averaging_operation_names ();
613 
619  AveragingOperation parse_averaging_operation_name (const std::string &s);
620 
626  template <int dim>
627  void average (const AveragingOperation operation,
628  const typename DoFHandler<dim>::active_cell_iterator &cell,
629  const Quadrature<dim> &quadrature_formula,
630  const Mapping<dim> &mapping,
631  MaterialModelOutputs<dim> &values_out);
632 
638  void average_property (const AveragingOperation operation,
639  const FullMatrix<double> &projection_matrix,
640  const FullMatrix<double> &expansion_matrix,
641  std::vector<double> &values_out);
642  }
643 
644 
663  template <int dim>
665  {
666  public:
668  {}
669 
674  virtual void
675  fill (const LinearAlgebra::BlockVector &solution,
676  const FEValuesBase<dim> &fe_values,
677  const Introspection<dim> &introspection) = 0;
678  };
679 
680 
707  template <int dim>
709  {
710  public:
715  virtual ~AdditionalMaterialOutputs() = default;
716 
717  virtual void average (const MaterialAveraging::AveragingOperation /*operation*/,
718  const FullMatrix<double> &/*projection_matrix*/,
719  const FullMatrix<double> &/*expansion_matrix*/)
720  {}
721  };
722 
723 
745  template <int dim>
747  {
748  public:
757  NamedAdditionalMaterialOutputs(const std::vector<std::string> &output_names);
758 
763 
768  const std::vector<std::string> &get_names() const;
769 
774  virtual std::vector<double> get_nth_output(const unsigned int idx) const = 0;
775 
776  virtual void average (const MaterialAveraging::AveragingOperation /*operation*/,
777  const FullMatrix<double> &/*projection_matrix*/,
778  const FullMatrix<double> &/*expansion_matrix*/)
779  {}
780 
781  private:
782  const std::vector<std::string> names;
783  };
784 
785 
791  template <int dim>
793  {
794  public:
795  SeismicAdditionalOutputs(const unsigned int n_points);
796 
797  virtual std::vector<double> get_nth_output(const unsigned int idx) const;
798 
804  std::vector<double> vs;
805 
811  std::vector<double> vp;
812  };
813 
814 
835  template <int dim>
837  {
838  public:
839  ReactionRateOutputs (const unsigned int n_points,
840  const unsigned int n_comp);
841 
842  virtual std::vector<double> get_nth_output(const unsigned int idx) const;
843 
851  std::vector<std::vector<double> > reaction_rates;
852  };
853 
854 
873  template <int dim>
875  {
876  public:
877  PrescribedFieldOutputs (const unsigned int n_points,
878  const unsigned int n_comp);
879 
880  virtual std::vector<double> get_nth_output(const unsigned int idx) const;
881 
889  std::vector<std::vector<double> > prescribed_field_outputs;
890  };
891 
892 
899  template <int dim>
901  {
902  public:
903  AdditionalMaterialOutputsStokesRHS(const unsigned int n_points)
904  : rhs_u(n_points), rhs_p(n_points), rhs_melt_pc(n_points)
905  {}
906 
908  {}
909 
910  virtual void average (const MaterialAveraging::AveragingOperation /*operation*/,
911  const FullMatrix<double> &/*projection_matrix*/,
912  const FullMatrix<double> &/*expansion_matrix*/)
913  {
914  // TODO: not implemented
915  }
916 
922  std::vector<Tensor<1,dim> > rhs_u;
923 
928  std::vector<double> rhs_p;
929 
934  std::vector<double> rhs_melt_pc;
935  };
936 
943  template <int dim>
945  {
946  public:
947  ElasticOutputs(const unsigned int n_points)
948  : elastic_force(n_points, numbers::signaling_nan<Tensor<2,dim> >() )
949  {}
950 
951  virtual ~ElasticOutputs()
952  {}
953 
954  virtual void average (const MaterialAveraging::AveragingOperation operation,
955  const FullMatrix<double> &/*projection_matrix*/,
956  const FullMatrix<double> &/*expansion_matrix*/)
957  {
959  return;
960  }
961 
967  std::vector<Tensor<2,dim> > elastic_force;
968  };
969 
970 
971 
987  template <int dim>
988  class Interface
989  {
990  public:
1005 
1010  virtual ~Interface();
1011 
1017  virtual
1018  void
1019  initialize ();
1020 
1025  virtual void update ();
1026 
1039  get_model_dependence () const;
1040 
1050  virtual bool is_compressible () const = 0;
1077  virtual double reference_viscosity () const = 0;
1087  virtual
1088  void evaluate (const MaterialModel::MaterialModelInputs<dim> &in,
1100  static
1101  void
1103 
1110  virtual
1111  void
1112  parse_parameters (ParameterHandler &prm);
1122  virtual
1123  void
1124  create_additional_named_outputs (MaterialModelOutputs &outputs) const;
1125 
1126 
1134  virtual
1135  void
1136  fill_additional_material_model_inputs(MaterialModel::MaterialModelInputs<dim> &input,
1137  const LinearAlgebra::BlockVector &solution,
1138  const FEValuesBase<dim> &fe_values,
1139  const Introspection<dim> &introspection) const;
1140 
1141  protected:
1158  };
1159 
1175  template <int dim>
1176  void
1177  register_material_model (const std::string &name,
1178  const std::string &description,
1179  void (*declare_parameters_function) (ParameterHandler &),
1180  Interface<dim> *(*factory_function) ());
1181 
1192  template <int dim>
1193  Interface<dim> *
1194  create_material_model (const std::string &model_name);
1195 
1196 
1208  template <int dim>
1209  Interface<dim> *
1211 
1212 
1219  template <int dim>
1220  std::string
1222 
1223 
1229  template <int dim>
1230  void
1232 
1233 
1234 
1244  template <int dim>
1245  void
1246  write_plugin_graph (std::ostream &output_stream);
1247 
1248 
1249 
1250 // --------------------- template function definitions ----------------------------------
1251 
1252  template <int dim>
1253  template <class AdditionalInputType>
1255  {
1256  for (unsigned int i=0; i<additional_inputs.size(); ++i)
1257  {
1258  AdditionalInputType *result = dynamic_cast<AdditionalInputType *> (additional_inputs[i].get());
1259  if (result)
1260  return result;
1261  }
1262  return nullptr;
1263  }
1264 
1265 
1266  template <int dim>
1267  template <class AdditionalInputType>
1268  const AdditionalInputType *MaterialModelInputs<dim>::get_additional_input() const
1269  {
1270  for (unsigned int i=0; i<additional_inputs.size(); ++i)
1271  {
1272  const AdditionalInputType *result = dynamic_cast<const AdditionalInputType *> (additional_inputs[i].get());
1273  if (result)
1274  return result;
1275  }
1276  return nullptr;
1277  }
1278 
1279 
1280  template <int dim>
1281  template <class AdditionalOutputType>
1283  {
1284  for (unsigned int i=0; i<additional_outputs.size(); ++i)
1285  {
1286  AdditionalOutputType *result = dynamic_cast<AdditionalOutputType *> (additional_outputs[i].get());
1287  if (result)
1288  return result;
1289  }
1290  return nullptr;
1291  }
1292 
1293 
1294  template <int dim>
1295  template <class AdditionalOutputType>
1296  const AdditionalOutputType *MaterialModelOutputs<dim>::get_additional_output() const
1297  {
1298  for (unsigned int i=0; i<additional_outputs.size(); ++i)
1299  {
1300  const AdditionalOutputType *result = dynamic_cast<const AdditionalOutputType *> (additional_outputs[i].get());
1301  if (result)
1302  return result;
1303  }
1304  return nullptr;
1305  }
1306 
1307 
1315 #define ASPECT_REGISTER_MATERIAL_MODEL(classname,name,description) \
1316  template class classname<2>; \
1317  template class classname<3>; \
1318  namespace ASPECT_REGISTER_MATERIAL_MODEL_ ## classname \
1319  { \
1320  aspect::internal::Plugins::RegisterHelper<aspect::MaterialModel::Interface<2>,classname<2> > \
1321  dummy_ ## classname ## _2d (&aspect::MaterialModel::register_material_model<2>, \
1322  name, description); \
1323  aspect::internal::Plugins::RegisterHelper<aspect::MaterialModel::Interface<3>,classname<3> > \
1324  dummy_ ## classname ## _3d (&aspect::MaterialModel::register_material_model<3>, \
1325  name, description); \
1326  }
1327  }
1328 }
1329 
1330 
1331 #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:1254
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:717
std::vector< std::vector< double > > prescribed_field_outputs
Definition: interface.h:889
MaterialModel::MaterialModelOutputs< dim > MaterialModelOutputs
Definition: interface.h:1004
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:947
AdditionalOutputType * get_additional_output()
Definition: interface.h:1282
std::vector< std::vector< double > > reaction_rates
Definition: interface.h:851
void declare_parameters(ParameterHandler &prm)
bool identifies_single_variable(const Dependence dependence)
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:782
std::vector< std::vector< double > > reaction_terms
Definition: interface.h:516
MaterialModel::MaterialModelInputs< dim > MaterialModelInputs
Definition: interface.h:997
AveragingOperation parse_averaging_operation_name(const std::string &s)
DoFHandler< dim >::active_cell_iterator current_cell
Definition: interface.h:350
virtual void average(const MaterialAveraging::AveragingOperation, const FullMatrix< double > &, const FullMatrix< double > &)
Definition: interface.h:776
AdditionalMaterialOutputsStokesRHS(const unsigned int n_points)
Definition: interface.h:903
std::vector< std::unique_ptr< AdditionalMaterialInputs< dim > > > additional_inputs
Definition: interface.h:357
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:967
Definition: compat.h:37
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
virtual void average(const MaterialAveraging::AveragingOperation, const FullMatrix< double > &, const FullMatrix< double > &)
Definition: interface.h:910
std::vector< double > entropy_derivative_temperature
Definition: interface.h:478
NonlinearDependence::ModelDependence model_dependence
Definition: interface.h:1157
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
virtual void average(const MaterialAveraging::AveragingOperation operation, const FullMatrix< double > &, const FullMatrix< double > &)
Definition: interface.h:954
std::vector< std::unique_ptr< AdditionalMaterialOutputs< dim > > > additional_outputs
Definition: interface.h:523