Hermes2D  2.0
discrete_problem.h
1 
2 
3 
4 
5 
6 
7 
8 
9 
10 
11 
12 
13 
14 
15 
16 #ifndef __H2D_DISCRETE_PROBLEM_H
17 #define __H2D_DISCRETE_PROBLEM_H
18 
19 #include "hermes_common.h"
20 #include "adapt/adapt.h"
21 #include "graph.h"
22 #include "forms.h"
23 #include "weakform/weakform.h"
24 #include "function/function.h"
25 #include "neighbor.h"
26 #include "refinement_selectors/selector.h"
27 #include "exceptions.h"
28 #include "mixins2d.h"
29 
30 namespace Hermes
31 {
32  namespace Hermes2D
33  {
34  class PrecalcShapeset;
35 
38  class NeighborNode
39  {
40  private:
41  NeighborNode(NeighborNode* parent, unsigned int transformation);
42  ~NeighborNode();
43  void set_left_son(NeighborNode* left_son);
44  void set_right_son(NeighborNode* right_son);
45  void set_transformation(unsigned int transformation);
46  NeighborNode* get_left_son();
47  NeighborNode* get_right_son();
48  unsigned int get_transformation();
49  NeighborNode* parent;
50  NeighborNode* left_son;
51  NeighborNode* right_son;
52  unsigned int transformation;
53  template<typename Scalar> friend class DiscreteProblem;
54  template<typename Scalar> friend class KellyTypeAdapt;
55  };
56 
62  template<typename Scalar>
63  class HERMES_API DiscreteProblem : public DiscreteProblemInterface<Scalar>, public Hermes::Mixins::TimeMeasurable, public Hermes::Hermes2D::Mixins::SettableSpaces<Scalar>, public Hermes::Hermes2D::Mixins::StateQueryable
64  {
65  public:
67  DiscreteProblem(const WeakForm<Scalar>* wf, Hermes::vector<const Space<Scalar> *> spaces);
68 
70  DiscreteProblem(const WeakForm<Scalar>* wf, const Space<Scalar>* space);
71 
73  virtual bool isOkay() const;
74  virtual inline std::string getClassName() const { return "DiscreteProblem"; }
75 
77  void set_fvm();
78 
80  virtual void set_spaces(Hermes::vector<const Space<Scalar>*> spaces);
81  virtual void set_space(const Space<Scalar>* space);
82 
84  DiscreteProblem();
85 
87  virtual ~DiscreteProblem();
88 
90  inline void set_do_not_use_cache() { this->do_not_use_cache = true; }
91 
93  const WeakForm<Scalar>* get_weak_formulation() const;
94 
96  void set_weak_formulation(const WeakForm<Scalar>* wf);
97 
99  virtual Hermes::vector<const Space<Scalar>*> get_spaces() const;
100 
102  int get_num_dofs() const;
103 
105  bool is_matrix_free() const;
106 
108  virtual void set_time(double time);
109  virtual void set_time_step(double time_step);
110 
111  void delete_cache();
112 
122  void assemble(Scalar* coeff_vec, SparseMatrix<Scalar>* mat, Vector<Scalar>* rhs = NULL,
123  bool force_diagonal_blocks = false, Table* block_weights = NULL);
124 
127  void assemble(Scalar* coeff_vec, Vector<Scalar>* rhs = NULL,
128  bool force_diagonal_blocks = false, Table* block_weights = NULL);
129 
132  virtual void assemble(SparseMatrix<Scalar>* mat, Vector<Scalar>* rhs = NULL, bool force_diagonal_blocks = false,
133  Table* block_weights = NULL);
134 
138  void assemble(Vector<Scalar>* rhs = NULL, bool force_diagonal_blocks = false,
139  Table* block_weights = NULL);
140 
143  static int init_geometry_points(RefMap* reference_mapping, int order, Geom<double>*& geometry, double*& jacobian_x_weights);
144  static int init_surface_geometry_points(RefMap* reference_mapping, int& order, Traverse::State* current_state, Geom<double>*& geometry, double*& jacobian_x_weights);
145 
146  protected:
147  void init_assembling(Scalar* coeff_vec, PrecalcShapeset*** pss , PrecalcShapeset*** spss, RefMap*** refmaps, Solution<Scalar>*** u_ext, AsmList<Scalar>*** als, WeakForm<Scalar>** weakforms);
148 
149  void deinit_assembling(PrecalcShapeset*** pss , PrecalcShapeset*** spss, RefMap*** refmaps, Solution<Scalar>*** u_ext, AsmList<Scalar>*** als, WeakForm<Scalar>** weakforms);
150 
152  bool form_to_be_assembled(MatrixForm<Scalar>* form, Traverse::State* current_state);
153  bool form_to_be_assembled(MatrixFormVol<Scalar>* form, Traverse::State* current_state);
154  bool form_to_be_assembled(MatrixFormSurf<Scalar>* form, Traverse::State* current_state);
155  bool form_to_be_assembled(VectorForm<Scalar>* form, Traverse::State* current_state);
156  bool form_to_be_assembled(VectorFormVol<Scalar>* form, Traverse::State* current_state);
157  bool form_to_be_assembled(VectorFormSurf<Scalar>* form, Traverse::State* current_state);
158 
159  // Return scaling coefficient.
160  double block_scaling_coeff(MatrixForm<Scalar>* form) const;
161 
169  void create_sparse_structure();
170  void create_sparse_structure(SparseMatrix<Scalar>* mat, Vector<Scalar>* rhs = NULL);
171 
173  inline void set_RK(int original_spaces_count) { this->RungeKutta = true; RK_original_spaces_count = original_spaces_count; }
174 
177  bool state_needs_recalculation(AsmList<Scalar>** current_als, Traverse::State* current_state);
178 
180  void calculate_cache_records(PrecalcShapeset** current_pss, PrecalcShapeset** current_spss, RefMap** current_refmaps, Solution<Scalar>** current_u_ext, AsmList<Scalar>** current_als,
181  Traverse::State* current_state, AsmList<Scalar>** current_alsSurface, WeakForm<Scalar>* current_wf);
182 
184  void assemble_one_state(PrecalcShapeset** current_pss, PrecalcShapeset** current_spss, RefMap** current_refmaps, Solution<Scalar>** current_u_ext,
185  AsmList<Scalar>** current_als, Traverse::State* current_state, WeakForm<Scalar>* current_wf);
186 
188  void adjust_order_to_refmaps(Form<Scalar> *form, int& order, Hermes::Ord* o, RefMap** current_refmaps);
189 
191  int calc_order_matrix_form(MatrixForm<Scalar>* mfv, RefMap** current_refmaps, Solution<Scalar>** current_u_ext, Traverse::State* current_state);
192 
194  virtual void assemble_matrix_form(MatrixForm<Scalar>* form, int order, Func<double>** base_fns, Func<double>** test_fns, Func<Scalar>** ext, Func<Scalar>** u_ext,
195  AsmList<Scalar>* current_als_i, AsmList<Scalar>* current_als_j, Traverse::State* current_state, int n_quadrature_points, Geom<double>* geometry, double* jacobian_x_weights);
196 
198  int calc_order_vector_form(VectorForm<Scalar>* mfv, RefMap** current_refmaps, Solution<Scalar>** current_u_ext, Traverse::State* current_state);
199 
201  void assemble_vector_form(VectorForm<Scalar>* form, int order, Func<double>** test_fns, Func<Scalar>** ext, Func<Scalar>** u_ext,
202  AsmList<Scalar>* current_als, Traverse::State* current_state, int n_quadrature_points, Geom<double>* geometry, double* jacobian_x_weights);
203 
206  void init_ext_orders(Form<Scalar> *form, Func<Hermes::Ord>** oi, Func<Hermes::Ord>** oext, Solution<Scalar>** current_u_ext, Traverse::State* current_state);
209  void deinit_ext_orders(Form<Scalar> *form, Func<Hermes::Ord>** oi, Func<Hermes::Ord>** oext);
210 
212  void init();
213 
215  bool is_up_to_date() const;
216 
218  const WeakForm<Scalar>* wf;
219 
221  Hermes::vector<const Space<Scalar>*> spaces;
222 
223  Hermes::vector<unsigned int> spaces_first_dofs;
224 
226  int* sp_seq;
227 
229  int ndof;
230 
232  Geom<Hermes::Ord> geom_ord;
233 
235  static double fake_wt;
236 
239  bool is_fvm;
240 
242  bool is_linear;
243 
246  bool have_matrix;
247 
249  bool DG_matrix_forms_present;
250 
252  bool DG_vector_forms_present;
253 
255  bool RungeKutta;
256 
258  int RK_original_spaces_count;
259 
261  SparseMatrix<Scalar>* current_mat;
262  Vector<Scalar>* current_rhs;
263  bool current_force_diagonal_blocks;
264  Table* current_block_weights;
265 
267  class CacheRecordPerElement
268  {
269  public:
270  void clear();
271  int* asmlistIdx;
272  int asmlistCnt;
273  };
274 
275  class CacheRecordPerSubIdx
276  {
277  public:
278  CacheRecordPerSubIdx();
279  int nvert;
280  int order;
281  void clear();
282  int asmlistCnt;
283  Func<double>** fns;
284  Func<double>*** fnsSurface;
285  Geom<double>* geometry;
286  Geom<double>** geometrySurface;
287  double* jacobian_x_weights;
288  double** jacobian_x_weightsSurface;
289  int n_quadrature_points;
290  int* n_quadrature_pointsSurface;
291  int* orderSurface;
292  int* asmlistSurfaceCnt;
293  };
294 
295  std::map<uint64_t, CacheRecordPerSubIdx*>*** cache_records_sub_idx;
296  CacheRecordPerElement*** cache_records_element;
297  bool** cache_element_stored;
298  int cache_size;
299  bool do_not_use_cache;
300 
302  Hermes::Exceptions::Exception* caughtException;
303 
304 
306 
308  void assemble_one_DG_state(PrecalcShapeset** current_pss, PrecalcShapeset** current_spss, RefMap** current_refmaps, Solution<Scalar>** current_u_ext, AsmList<Scalar>** current_als,
309  Traverse::State* current_state, Hermes::vector<MatrixFormDG<Scalar>*> current_mfDG, Hermes::vector<VectorFormDG<Scalar>*> current_vfDG, Transformable** fn);
310 
312  void assemble_DG_one_neighbor(bool edge_processed, unsigned int neighbor_i,
313  PrecalcShapeset** current_pss, PrecalcShapeset** current_spss, RefMap** current_refmaps, Solution<Scalar>** current_u_ext, AsmList<Scalar>** current_als,
314  Traverse::State* current_state, Hermes::vector<MatrixFormDG<Scalar>*> current_mfDG, Hermes::vector<VectorFormDG<Scalar>*> current_vfDG, Transformable** fn,
315  std::map<unsigned int, PrecalcShapeset *> npss, std::map<unsigned int, PrecalcShapeset *> nspss, std::map<unsigned int, RefMap *> nrefmap,
316  LightArray<NeighborSearch<Scalar>*>& neighbor_searches, unsigned int min_dg_mesh_seq);
317 
319  void assemble_DG_matrix_forms(PrecalcShapeset** current_pss, PrecalcShapeset** current_spss, RefMap** current_refmaps, Solution<Scalar>** current_u_ext, AsmList<Scalar>** current_als,
320  Traverse::State* current_state, MatrixFormDG<Scalar>** current_mfDG, std::map<unsigned int, PrecalcShapeset*> npss,
321  std::map<unsigned int, PrecalcShapeset*> nspss, std::map<unsigned int, RefMap*> nrefmap, LightArray<NeighborSearch<Scalar>*>& neighbor_searches);
322 
324  void assemble_DG_vector_forms(PrecalcShapeset** current_spss, RefMap** current_refmaps, Solution<Scalar>** current_u_ext, AsmList<Scalar>** current_als,
325  Traverse::State* current_state, VectorFormDG<Scalar>** current_vfDG, std::map<unsigned int, PrecalcShapeset*> nspss,
326  std::map<unsigned int, RefMap*> nrefmap, LightArray<NeighborSearch<Scalar>*>& neighbor_searches);
327 
328  DiscontinuousFunc<Hermes::Ord>* init_ext_fn_ord(NeighborSearch<Scalar>* ns, MeshFunction<Scalar>* fu);
329 
331  int calc_order_dg_matrix_form(MatrixFormDG<Scalar>* mfDG, Hermes::vector<Solution<Scalar>*> u_ext,
332  PrecalcShapeset* fu, PrecalcShapeset* fv, RefMap* ru, SurfPos* surf_pos,
333  bool neighbor_supp_u, bool neighbor_supp_v, LightArray<NeighborSearch<Scalar>*>& neighbor_searches, int neighbor_index_u);
334 
336  Scalar eval_dg_form(MatrixFormDG<Scalar>* mfDG, Hermes::vector<Solution<Scalar>*> u_ext,
337  PrecalcShapeset* fu, PrecalcShapeset* fv, RefMap* ru_central, RefMap* ru_actual, RefMap* rv,
338  bool neighbor_supp_u, bool neighbor_supp_v,
339  SurfPos* surf_pos, LightArray<NeighborSearch<Scalar>*>& neighbor_searches, int neighbor_index_u, int neighbor_index_v);
340 
342  int calc_order_dg_vector_form(VectorFormDG<Scalar>* vfDG, Hermes::vector<Solution<Scalar>*> u_ext,
343  PrecalcShapeset* fv, RefMap* ru, SurfPos* surf_pos,
344  LightArray<NeighborSearch<Scalar>*>& neighbor_searches, int neighbor_index_v);
345 
347  Scalar eval_dg_form(VectorFormDG<Scalar>* vfDG, Hermes::vector<Solution<Scalar>*> u_ext,
348  PrecalcShapeset* fv, RefMap* rv,
349  SurfPos* surf_pos, LightArray<NeighborSearch<Scalar>*>& neighbor_searches, int neighbor_index_v);
350 
352  Func<Hermes::Ord>** init_ext_fns_ord(Hermes::vector<MeshFunction<Scalar>*> &ext,
353  LightArray<NeighborSearch<Scalar>*>& neighbor_searches);
354 
356  Func<Scalar>** init_ext_fns(Hermes::vector<MeshFunction<Scalar>*> &ext,
357  LightArray<NeighborSearch<Scalar>*>& neighbor_searches,
358  int order, unsigned int min_dg_mesh_seq);
359 
361  bool init_neighbors(LightArray<NeighborSearch<Scalar>*>& neighbor_searches, Traverse::State* current_state, unsigned int min_dg_mesh_seq);
362 
364  void build_multimesh_tree(NeighborNode* root, LightArray<NeighborSearch<Scalar>*>& neighbor_searches);
365 
367  void insert_into_multimesh_tree(NeighborNode* node, unsigned int* transformations, unsigned int transformation_count);
368 
370  Hermes::vector<Hermes::vector<unsigned int>*> get_multimesh_neighbors_transformations(NeighborNode* multimesh_tree);
371 
373  void traverse_multimesh_tree(NeighborNode* node, Hermes::vector<Hermes::vector<unsigned int>*>& running_transformations);
374 
376  void update_neighbor_search(NeighborSearch<Scalar>* ns, NeighborNode* multimesh_tree);
377 
380  NeighborNode* find_node(unsigned int* transformations, unsigned int transformation_count, NeighborNode* node);
381 
384  int update_ns_subtree(NeighborSearch<Scalar>* ns, NeighborNode* node, unsigned int ith_neighbor);
385 
387  void traverse_multimesh_subtree(NeighborNode* node, Hermes::vector<Hermes::vector<unsigned int>*>& running_central_transformations,
388  Hermes::vector<Hermes::vector<unsigned int>*>& running_neighbor_transformations, const typename NeighborSearch<Scalar>::NeighborEdgeInfo& edge_info, const int& active_edge, const int& mode);
389 
390  template<typename T> friend class KellyTypeAdapt;
391  template<typename T> friend class LinearSolver;
392  template<typename T> friend class NewtonSolver;
393  template<typename T> friend class PicardSolver;
394  template<typename T> friend class RungeKutta;
395  };
396  }
397 }
398 #endif