16 #include "discrete_problem/dg/discrete_problem_dg_assembler.h"
17 #include "discrete_problem/discrete_problem_helpers.h"
18 #include "discrete_problem/discrete_problem.h"
19 #include "function/exact_solution.h"
20 #include "mesh/traverse.h"
21 #include "space/space.h"
22 #include "function/solution.h"
31 template<
typename Scalar>
34 this->init(to_set, dirichlet_lift_accordingly);
35 this->set_spaces(spaces);
36 this->set_weak_formulation(wf_);
39 template<
typename Scalar>
42 this->init(to_set, dirichlet_lift_accordingly);
43 this->set_space(space);
44 this->set_weak_formulation(wf_);
47 template<
typename Scalar>
50 init(to_set, dirichlet_lift_accordingly);
53 template<
typename Scalar>
56 this->reassembled_states_reuse_linear_system =
nullptr;
58 this->spaces_size = this->spaces.size();
60 this->nonlinear = !to_set;
61 if (dirichlet_lift_accordingly)
62 this->add_dirichlet_lift = !this->nonlinear;
64 this->add_dirichlet_lift = this->nonlinear;
66 if (this->add_dirichlet_lift)
67 this->dirichlet_lift_rhs = create_vector<Scalar>(
false);
69 this->dirichlet_lift_rhs =
nullptr;
74 for (
int i = 0; i < this->num_threads_used; i++)
78 template<
typename Scalar>
81 for (
int i = 0; i < this->num_threads_used; i++)
82 delete this->threadAssembler[i];
83 delete[] this->threadAssembler;
85 if (this->dirichlet_lift_rhs)
86 delete this->dirichlet_lift_rhs;
89 template<
typename Scalar>
95 if (this->spaces_size == 0)
99 for (
unsigned short space_i = 0; space_i < this->spaces_size; space_i++)
100 this->spaces[space_i]->check();
102 for (
unsigned short space_i = 0; space_i < this->spaces_size; space_i++)
103 if (!this->spaces[space_i]->is_up_to_date())
104 throw Exceptions::Exception(
"Space is out of date, if you manually refine it, you have to call assign_dofs().");
109 template<
typename Scalar>
112 Loggable::set_verbose_output(to_set);
113 this->selectiveAssembler.set_verbose_output(to_set);
116 template<
typename Scalar>
120 this->wf->set_current_time(time);
123 template<
typename Scalar>
126 this->wf->set_current_time_step(time_step);
129 template<
typename Scalar>
135 template<
typename Scalar>
140 this->selectiveAssembler.set_RK(original_spaces_count, force_diagonal_blocks_, block_weights_);
142 for (
int i = 0; i < this->num_threads_used; i++)
143 this->threadAssembler[i]->set_RK(original_spaces_count, force_diagonal_blocks_, block_weights_);
146 template<
typename Scalar>
149 this->selectiveAssembler.matrix_structure_reusable =
false;
152 template<
typename Scalar>
157 this->selectiveAssembler.set_weak_formulation(wf);
158 this->selectiveAssembler.matrix_structure_reusable =
false;
161 template<
typename Scalar>
166 for (
int i = 0; i < this->num_threads_used; i++)
167 this->threadAssembler[i]->set_matrix(mat);
169 if (mat && this->current_mat != mat)
171 this->invalidate_matrix();
178 template<
typename Scalar>
183 for (
int i = 0; i < this->num_threads_used; i++)
185 this->threadAssembler[i]->set_rhs(rhs);
186 this->threadAssembler[i]->dirichlet_lift_rhs = this->dirichlet_lift_rhs;
189 if (rhs && this->current_rhs != rhs)
191 this->invalidate_matrix();
198 template<
typename Scalar>
201 if (this->spaces_size > 0)
202 Helpers::check_length(spacesToSet, this->spaces_size);
204 for (
unsigned int i = 0; i < spacesToSet.size(); i++)
207 throw Exceptions::NullException(0, i);
208 spacesToSet[i]->check();
211 this->spaces_size = spacesToSet.size();
212 this->spaces = spacesToSet;
214 this->selectiveAssembler.set_spaces(spacesToSet);
216 for (
int i = 0; i < this->num_threads_used; i++)
217 this->threadAssembler[i]->init_spaces(spaces);
220 template<
typename Scalar>
223 std::vector<SpaceSharedPtr<Scalar> > spaces;
224 spaces.push_back(space);
225 this->set_spaces(spaces);
228 template<
typename Scalar>
231 Scalar* coeff_vec =
nullptr;
232 return assemble(coeff_vec, mat, rhs);
235 template<
typename Scalar>
238 return assemble(coeff_vec,
nullptr, rhs);
241 template<
typename Scalar>
244 Scalar* coeff_vec =
nullptr;
245 return assemble(coeff_vec,
nullptr, rhs);
248 template<
typename Scalar>
252 for (
unsigned int space_i = 0; space_i < spaces.size(); space_i++)
253 meshes.push_back(spaces[space_i]->get_mesh());
254 for (
unsigned int ext_i = 0; ext_i < this->wf->ext.size(); ext_i++)
255 meshes.push_back(this->wf->ext[ext_i]->get_mesh());
256 for (
unsigned int form_i = 0; form_i < this->wf->get_forms().size(); form_i++)
257 for (
unsigned int ext_i = 0; ext_i < this->wf->get_forms()[form_i]->ext.size(); ext_i++)
258 if (this->wf->get_forms()[form_i]->ext[ext_i])
259 meshes.push_back(this->wf->get_forms()[form_i]->ext[ext_i]->get_mesh());
263 for (
unsigned int space_i = 0; space_i < spaces.size(); space_i++)
264 meshes.push_back(spaces[space_i]->get_mesh());
271 for (
unsigned char i = 0; i < this->num_threads_used; i++)
272 this->threadAssembler[i]->set_weak_formulation(this->wf);
278 this->exceptionMessageCaughtInParallelBlock.clear();
281 if (this->add_dirichlet_lift)
284 this->dirichlet_lift_rhs->alloc(ndof);
288 template<
typename Scalar>
296 bool result = this->set_matrix(mat) && this->set_rhs(rhs);
299 unsigned int num_states;
301 std::vector<MeshSharedPtr> meshes;
302 this->init_assembling(states, num_states, meshes);
304 this->info(
"\tDiscreteProblem: Initialization: %s.", this->last_str().c_str());
309 if (this->selectiveAssembler.prepare_sparse_structure(this->current_mat, this->current_rhs, this->spaces, states, num_states))
312 this->info(
"\tDiscreteProblem: Prepare sparse structure: %s.", this->last_str().c_str());
315 if (this->current_mat && this->reassembled_states_reuse_linear_system)
316 this->reassembled_states_reuse_linear_system(states, num_states, this->current_mat, this->current_rhs, this->dirichlet_lift_rhs, coeff_vec);
319 if (this->nonlinear && coeff_vec)
323 for (
int i = 0; i < this->spaces_size; i++)
327 first_dof += spaces[i]->get_num_dofs();
334 bool is_DG = this->wf->is_DG();
336 #pragma omp parallel num_threads(this->num_threads_used)
338 int thread_number = omp_get_thread_num();
339 int start = (num_states / this->num_threads_used) * thread_number;
340 int end = (num_states / this->num_threads_used) * (thread_number + 1);
341 if (thread_number == this->num_threads_used - 1)
346 this->threadAssembler[thread_number]->init_assembling(u_ext_sln, spaces, this->add_dirichlet_lift);
352 for (
int state_i = start; state_i < end; state_i++)
355 if (!this->exceptionMessageCaughtInParallelBlock.empty())
376 this->threadAssembler[thread_number]->deinit_assembling();
378 catch (Hermes::Exceptions::Exception& e)
380 #pragma omp critical (exceptionMessageCaughtInParallelBlock)
381 this->exceptionMessageCaughtInParallelBlock = e.info();
385 #pragma omp critical (exceptionMessageCaughtInParallelBlock)
386 this->exceptionMessageCaughtInParallelBlock = e.what();
391 if (this->nonlinear && coeff_vec)
393 for (
int i = 0; i < this->spaces_size; i++)
402 this->deinit_assembling(states, num_states);
405 if (this->current_mat)
406 this->current_mat->finish();
407 if (this->current_rhs)
408 this->current_rhs->finish();
410 if (!this->exceptionMessageCaughtInParallelBlock.empty())
411 throw Hermes::Exceptions::Exception(this->exceptionMessageCaughtInParallelBlock.c_str());
414 for (
unsigned int space_i = 0; space_i < spaces.size(); space_i++)
416 for_all_active_elements(e, spaces[space_i]->get_mesh())
418 spaces[space_i]->edata[e->
id].changed_in_last_adaptation =
false;
424 this->info(
"\tDiscreteProblem: De-initialization: %s.", this->last_str().c_str());
429 template<
typename Scalar>
432 for (
unsigned int i = 0; i < num_states; i++)
434 free_with_check(states);
437 if (this->add_dirichlet_lift && this->current_rhs)
438 this->current_rhs->add_vector(this->dirichlet_lift_rhs);
441 template class HERMES_API DiscreteProblem < double > ;
442 template class HERMES_API DiscreteProblem < std::complex<double> > ;
Stores one element of a mesh.
bool visited
true if the element has been visited during assembling
::xsd::cxx::tree::exception< char > exception
Root of the C++/Tree exception hierarchy.
Used to pass the instances of Space around.
::xsd::cxx::tree::time< char, simple_type > time
C++ type corresponding to the time XML Schema built-in type.
void deinit_assembling_one_state()
Deinitialize assembling for a state.
void assemble_one_state()
Assemble DG forms.
This class is a one-thread (non-DG) assembly worker.
Represents a finite element space over a domain.
void init_assembling_one_state(Traverse::State *current_state_)
Initialize assembling for a state.
State ** get_states(std::vector< MeshSharedPtr > meshes, unsigned int &states_count)
Represents the solution of a PDE.