Hermes2D  2.0
ogprojection.cpp
1 // This file is part of Hermes2D.
2 //
3 // Hermes2D is free software: you can redistribute it and/or modify
4 // it under the terms of the GNU General Public License as published by
5 // the Free Software Foundation, either version 2 of the License, or
6 // (at your option) any later version.
7 //
8 // Hermes2D is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.See the
11 // GNU General Public License for more details.
12 //
13 // You should have received a copy of the GNU General Public License
14 // along with Hermes2D. If not, see <http://www.gnu.org/licenses/>.
15 
16 #include "projections/ogprojection.h"
17 #include "space.h"
18 #include "linear_solver.h"
19 
20 namespace Hermes
21 {
22  namespace Hermes2D
23  {
24  template<typename Scalar>
25  OGProjection<Scalar>::OGProjection() : Hermes::Mixins::Loggable(), ndof(0)
26  {
27  }
28 
29  template<typename Scalar>
30  void OGProjection<Scalar>::project_internal(const Space<Scalar>* space, WeakForm<Scalar>* wf,
31  Scalar* target_vec)
32  {
33  // Sanity check.
34  if(space == NULL)
35  throw Hermes::Exceptions::Exception("this->space == NULL in project_internal().");
36 
37  // Initialize DiscreteProblem.
38  DiscreteProblemLinear<Scalar> dp(wf, space);
39  dp.set_do_not_use_cache();
40 
41  // Initialize linear solver.
42  Hermes::Hermes2D::LinearSolver<Scalar> linear_solver(&dp);
43  linear_solver.set_verbose_output(this->get_verbose_output());
44 
45  // Perform Newton iteration.
46  linear_solver.solve();
47 
48  if(target_vec != NULL)
49  for (int i = 0; i < space->get_num_dofs(); i++)
50  target_vec[i] = linear_solver.get_sln_vector()[i];
51  }
52 
53  template<typename Scalar>
54  void OGProjection<Scalar>::project_global(const Space<Scalar>* space,
55  MatrixFormVol<Scalar>* custom_projection_jacobian,
56  VectorFormVol<Scalar>* custom_projection_residual,
57  Scalar* target_vec)
58  {
59  // Define projection weak form.
60  WeakForm<Scalar>* proj_wf = new WeakForm<Scalar>(1);
61  proj_wf->add_matrix_form(custom_projection_jacobian);
62  proj_wf->add_vector_form(custom_projection_residual);
63 
64  // Call the main function.
65  project_internal(space, proj_wf, target_vec);
66 
67  // Clean up.
68  delete proj_wf;
69  }
70 
71  template<typename Scalar>
72  void OGProjection<Scalar>::project_global(const Space<Scalar>* space,
73  MeshFunction<Scalar>* source_meshfn, Scalar* target_vec,
74  ProjNormType proj_norm)
75  {
76  // Sanity checks.
77  if(target_vec == NULL)
78  throw Exceptions::NullException(3);
79 
80  // If projection norm is not provided, set it
81  // to match the type of the space.
82  ProjNormType norm = HERMES_UNSET_NORM;
83  if(proj_norm == HERMES_UNSET_NORM)
84  {
85  SpaceType space_type = space->get_type();
86  switch (space_type)
87  {
88  case HERMES_H1_SPACE: norm = HERMES_H1_NORM; break;
89  case HERMES_HCURL_SPACE: norm = HERMES_HCURL_NORM; break;
90  case HERMES_HDIV_SPACE: norm = HERMES_HDIV_NORM; break;
91  case HERMES_L2_SPACE: norm = HERMES_L2_NORM; break;
92  default: throw Hermes::Exceptions::Exception("Unknown space type in OGProjection<Scalar>::project_global().");
93  }
94  }
95  else norm = proj_norm;
96 
97  // Define temporary projection weak form.
98  WeakForm<Scalar>* proj_wf = new WeakForm<Scalar>(1);
99  proj_wf->warned_nonOverride = true;
100  proj_wf->set_ext(source_meshfn);
101  // Add Jacobian.
102  proj_wf->add_matrix_form(new ProjectionMatrixFormVol(0, 0, norm));
103  // Add Residual.
104  proj_wf->add_vector_form(new ProjectionVectorFormVol(0, norm));
105 
106  // Call main function.
107  project_internal(space, proj_wf, target_vec);
108 
109  // Clean up.
110  delete proj_wf;
111  }
112 
113  template<typename Scalar>
114  void OGProjection<Scalar>::project_global(const Space<Scalar>* space,
115  Solution<Scalar>* source_sln, Solution<Scalar>* target_sln,
116  ProjNormType proj_norm)
117  {
118  if(proj_norm == HERMES_UNSET_NORM)
119  {
120  SpaceType space_type = space->get_type();
121  switch (space_type)
122  {
123  case HERMES_H1_SPACE: proj_norm = HERMES_H1_NORM; break;
124  case HERMES_HCURL_SPACE: proj_norm = HERMES_HCURL_NORM; break;
125  case HERMES_HDIV_SPACE: proj_norm = HERMES_HDIV_NORM; break;
126  case HERMES_L2_SPACE: proj_norm = HERMES_L2_NORM; break;
127  default: throw Hermes::Exceptions::Exception("Unknown space type in OGProjection<Scalar>::project_global().");
128  }
129  }
130 
131  // Calculate the coefficient vector.
132  int ndof = space->get_num_dofs();
133  Scalar* target_vec = new Scalar[ndof];
134  project_global(space, source_sln, target_vec, proj_norm);
135 
136  // Translate coefficient vector into a Solution.
137  Solution<Scalar>::vector_to_solution(target_vec, space, target_sln);
138 
139  // Clean up.
140  delete [] target_vec;
141  }
142 
143  template<typename Scalar>
144  void OGProjection<Scalar>::project_global(Hermes::vector<const Space<Scalar>*> spaces,
145  Hermes::vector<MeshFunction<Scalar>*> source_meshfns,
146  Scalar* target_vec, Hermes::vector<ProjNormType> proj_norms)
147  {
148  int n = spaces.size();
149 
150  // Sanity checks.
151  if(n != source_meshfns.size())
152  throw Exceptions::LengthException(1, 2, n, source_meshfns.size());
153  if(target_vec == NULL)
154  throw Exceptions::NullException(3);
155  if(!proj_norms.empty() && n != proj_norms.size())
156  throw Exceptions::LengthException(1, 5, n, proj_norms.size());
157 
158  int start_index = 0;
159  for (int i = 0; i < n; i++)
160  {
161  if(i == 0)
162  this->info("Projection: %d-th space", i);
163  if(i == 1)
164  this->info("Projection: %d-st space", i);
165  if(i == 2)
166  this->info("Projection: %d-nd space", i);
167  if(i == 3)
168  this->info("Projection: %d-rd space", i);
169  if(i > 3)
170  this->info("Projection: %d-th space", i);
171  if(proj_norms.empty())
172  project_global(spaces[i], source_meshfns[i], target_vec + start_index, HERMES_UNSET_NORM);
173  else
174  project_global(spaces[i], source_meshfns[i], target_vec + start_index, proj_norms[i]);
175  spaces[i]->assign_dofs(start_index);
176  start_index += spaces[i]->get_num_dofs();
177  }
178  }
179 
180  template<typename Scalar>
181  void OGProjection<Scalar>::project_global(Hermes::vector<const Space<Scalar>*> spaces, Hermes::vector<Solution<Scalar>*> source_slns,
182  Scalar* target_vec, Hermes::vector<ProjNormType> proj_norms)
183  {
184  int n = spaces.size();
185 
186  // Sanity checks.
187  if(n != source_slns.size())
188  throw Exceptions::LengthException(1, 2, n, source_slns.size());
189  if(target_vec == NULL)
190  throw Exceptions::NullException(3);
191  if(!proj_norms.empty() && n != proj_norms.size())
192  throw Exceptions::LengthException(1, 5, n, proj_norms.size());
193 
194  int start_index = 0;
195  for (int i = 0; i < n; i++)
196  {
197  if(i == 0)
198  this->info("Projection: %d-th space", i);
199  if(i == 1)
200  this->info("Projection: %d-st space", i);
201  if(i == 2)
202  this->info("Projection: %d-nd space", i);
203  if(i == 3)
204  this->info("Projection: %d-rd space", i);
205  if(i > 3)
206  this->info("Projection: %d-th space", i);
207  if(proj_norms.empty())
208  project_global(spaces[i], source_slns[i], target_vec + start_index, HERMES_UNSET_NORM);
209  else
210  project_global(spaces[i], source_slns[i], target_vec + start_index, proj_norms[i]);
211  start_index += spaces[i]->get_num_dofs();
212  }
213  }
214 
215  template<typename Scalar>
216  void OGProjection<Scalar>::project_global(Hermes::vector<const Space<Scalar>*> spaces, Hermes::vector<Solution<Scalar>*> source_slns,
217  Hermes::vector<Solution<Scalar>*> target_slns, Hermes::vector<ProjNormType> proj_norms, bool delete_old_meshes)
218  {
219  int n = spaces.size();
220 
221  // Sanity checks.
222  if(n != source_slns.size())
223  throw Exceptions::LengthException(1, 2, n, source_slns.size());
224  if(n != target_slns.size())
225  throw Exceptions::LengthException(1, 2, n, target_slns.size());
226  if(!proj_norms.empty() && n != proj_norms.size())
227  throw Exceptions::LengthException(1, 5, n, proj_norms.size());
228 
229  int start_index = 0;
230  for (int i = 0; i < n; i++)
231  {
232  if(i == 0)
233  this->info("Projection: %d-th space", i);
234  if(i == 1)
235  this->info("Projection: %d-st space", i);
236  if(i == 2)
237  this->info("Projection: %d-nd space", i);
238  if(i == 3)
239  this->info("Projection: %d-rd space", i);
240  if(i > 3)
241  this->info("Projection: %d-th space", i);
242  if(proj_norms.empty())
243  project_global(spaces[i], source_slns[i], target_slns[i], HERMES_UNSET_NORM);
244  else
245  project_global(spaces[i], source_slns[i], target_slns[i], proj_norms[i]);
246  start_index += spaces[i]->get_num_dofs();
247  }
248  }
249 
250  template class HERMES_API OGProjection<double>;
251  template class HERMES_API OGProjection<std::complex<double> >;
252  }
253 }