mesh_function.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // C++ includes
21 
22 
23 // Local Includes
24 #include "libmesh/mesh_function.h"
25 #include "libmesh/dense_vector.h"
27 #include "libmesh/numeric_vector.h"
28 #include "libmesh/dof_map.h"
30 #include "libmesh/fe_base.h"
31 #include "libmesh/fe_interface.h"
33 #include "libmesh/mesh_base.h"
34 #include "libmesh/point.h"
35 
36 namespace libMesh
37 {
38 
39 
40 //------------------------------------------------------------------
41 // MeshFunction methods
43  const NumericVector<Number>& vec,
44  const DofMap& dof_map,
45  const std::vector<unsigned int>& vars,
46  const FunctionBase<Number>* master) :
47  FunctionBase<Number> (master),
48  ParallelObject (eqn_systems),
49  _eqn_systems (eqn_systems),
50  _vector (vec),
51  _dof_map (dof_map),
52  _system_vars (vars),
53  _point_locator (NULL),
54  _out_of_mesh_mode (false),
55  _out_of_mesh_value ()
56 {
57 }
58 
59 
60 
62  const NumericVector<Number>& vec,
63  const DofMap& dof_map,
64  const unsigned int var,
65  const FunctionBase<Number>* master) :
66  FunctionBase<Number> (master),
67  ParallelObject (eqn_systems),
68  _eqn_systems (eqn_systems),
69  _vector (vec),
70  _dof_map (dof_map),
71  _system_vars (1,var),
72  _point_locator (NULL),
73  _out_of_mesh_mode (false),
74  _out_of_mesh_value ()
75 {
76 // std::vector<unsigned int> buf (1);
77 // buf[0] = var;
78 // _system_vars (buf);
79 }
80 
81 
82 
83 
84 
85 
86 
88 {
89  // only delete the point locator when we are the master
90  if ((this->_point_locator != NULL) && (this->_master == NULL))
91  delete this->_point_locator;
92 }
93 
94 
95 
96 
97 void MeshFunction::init (const Trees::BuildType /*point_locator_build_type*/)
98 {
99  // are indices of the desired variable(s) provided?
100  libmesh_assert_greater (this->_system_vars.size(), 0);
101 
102  // Don't do twice...
103  if (this->_initialized)
104  {
106  return;
107  }
108 
109  /*
110  * set up the PointLocator: either someone else
111  * is the master (go and get the address of his
112  * point locator) or this object is the master
113  * (build the point locator on our own).
114  */
115  if (this->_master != NULL)
116  {
117  // we aren't the master
118  const MeshFunction* master =
119  libmesh_cast_ptr<const MeshFunction*>(this->_master);
120 
121  if (master->_point_locator == NULL)
122  {
123  libMesh::err << "ERROR: When the master-servant concept is used,"
124  << std::endl
125  << " the master has to be initialized first!"
126  << std::endl;
127  libmesh_error();
128  }
129  else
130  {
131  this->_point_locator = master->_point_locator;
132  }
133  }
134  else
135  {
136  // we are the master: build the point locator
137 
138  // constant reference to the other mesh
139  const MeshBase& mesh = this->_eqn_systems.get_mesh();
140 
141  // build the point locator. Only \p TREE version available
142  //AutoPtr<PointLocatorBase> ap (PointLocatorBase::build (TREE, mesh));
143  //this->_point_locator = ap.release();
144  // this->_point_locator = new PointLocatorTree (mesh, point_locator_build_type);
145  this->_point_locator = mesh.sub_point_locator().release();
146 
147  // Point locator no longer needs to be initialized.
148  // this->_point_locator->init();
149  }
150 
151 
152  // ready for use
153  this->_initialized = true;
154 }
155 
156 
157 void
159 {
160  // only delete the point locator when we are the master
161  if ((this->_point_locator != NULL) && (this->_master == NULL))
162  {
163  delete this->_point_locator;
164  this->_point_locator = NULL;
165  }
166  this->_initialized = false;
167 }
168 
169 
170 
172 {
174  (new MeshFunction
176 }
177 
178 
179 
181  const Real time)
182 {
183  libmesh_assert (this->initialized());
184 
185  DenseVector<Number> buf (1);
186  this->operator() (p, time, buf);
187  return buf(0);
188 }
189 
190 
191 
193  const Real time)
194 {
195  libmesh_assert (this->initialized());
196 
197  std::vector<Gradient> buf (1);
198  this->gradient(p, time, buf);
199  return buf[0];
200 }
201 
202 
203 
204 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
206  const Real time)
207 {
208  libmesh_assert (this->initialized());
209 
210  std::vector<Tensor> buf (1);
211  this->hessian(p, time, buf);
212  return buf[0];
213 }
214 #endif
215 
216 
217 
219  const Real,
220  DenseVector<Number>& output)
221 {
222  libmesh_assert (this->initialized());
223 
224  /* Ensure that in the case of a master mesh function, the
225  out-of-mesh mode is enabled either for both or for none. This is
226  important because the out-of-mesh mode is also communicated to
227  the point locator. Since this is time consuming, enable it only
228  in debug mode. */
229 #ifdef DEBUG
230  if (this->_master != NULL)
231  {
232  const MeshFunction* master =
233  libmesh_cast_ptr<const MeshFunction*>(this->_master);
235  {
236  libMesh::err << "ERROR: If you use out-of-mesh-mode in connection with master mesh functions, you must enable out-of-mesh mode for both the master and the slave mesh function." << std::endl;
237  libmesh_error();
238  }
239  }
240 #endif
241 
242  // locate the point in the other mesh
243  const Elem* element = this->_point_locator->operator()(p);
244 
245  // If we have an element, but it's not a local element, then we
246  // either need to have a serialized vector or we need to find a
247  // local element sharing the same point.
248  if (element &&
249  (element->processor_id() != this->processor_id()) &&
250  _vector.type() != SERIAL)
251  {
252  // look for a local element containing the point
253  std::set<const Elem*> point_neighbors;
254  element->find_point_neighbors(p, point_neighbors);
255  element = NULL;
256  std::set<const Elem*>::const_iterator it = point_neighbors.begin();
257  const std::set<const Elem*>::const_iterator end = point_neighbors.end();
258  for (; it != end; ++it)
259  {
260  const Elem* elem = *it;
261  if (elem->processor_id() == this->processor_id())
262  {
263  element = elem;
264  break;
265  }
266  }
267  }
268 
269  if (!element)
270  {
271  output = _out_of_mesh_value;
272  }
273  else
274  {
275  // resize the output vector to the number of output values
276  // that the user told us
277  output.resize (libmesh_cast_int<unsigned int>
278  (this->_system_vars.size()));
279 
280 
281  {
282  const unsigned int dim = this->_eqn_systems.get_mesh().mesh_dimension();
283 
284 
285  /*
286  * Get local coordinates to feed these into compute_data().
287  * Note that the fe_type can safely be used from the 0-variable,
288  * since the inverse mapping is the same for all FEFamilies
289  */
290  const Point mapped_point (FEInterface::inverse_map (dim,
291  this->_dof_map.variable_type(0),
292  element,
293  p));
294 
295 
296  // loop over all vars
297  for (unsigned int index=0; index < this->_system_vars.size(); index++)
298  {
299  /*
300  * the data for this variable
301  */
302  const unsigned int var = _system_vars[index];
303  const FEType& fe_type = this->_dof_map.variable_type(var);
304 
309  {
310  FEComputeData data (this->_eqn_systems, mapped_point);
311 
312  FEInterface::compute_data (dim, fe_type, element, data);
313 
314  // where the solution values for the var-th variable are stored
315  std::vector<dof_id_type> dof_indices;
316  this->_dof_map.dof_indices (element, dof_indices, var);
317 
318  // interpolate the solution
319  {
320  Number value = 0.;
321 
322  for (unsigned int i=0; i<dof_indices.size(); i++)
323  value += this->_vector(dof_indices[i]) * data.shape[i];
324 
325  output(index) = value;
326  }
327 
328  }
329 
330  // next variable
331  }
332  }
333  }
334 
335  // all done
336  return;
337 }
338 
339 
340 
342  const Real,
343  std::vector<Gradient>& output)
344 {
345  libmesh_assert (this->initialized());
346 
347  /* Ensure that in the case of a master mesh function, the
348  out-of-mesh mode is enabled either for both or for none. This is
349  important because the out-of-mesh mode is also communicated to
350  the point locator. Since this is time consuming, enable it only
351  in debug mode. */
352 #ifdef DEBUG
353  if (this->_master != NULL)
354  {
355  const MeshFunction* master =
356  libmesh_cast_ptr<const MeshFunction*>(this->_master);
358  {
359  libMesh::err << "ERROR: If you use out-of-mesh-mode in connection with master mesh functions, you must enable out-of-mesh mode for both the master and the slave mesh function." << std::endl;
360  libmesh_error();
361  }
362  }
363 #endif
364 
365  // locate the point in the other mesh
366  const Elem* element = this->_point_locator->operator()(p);
367 
368  // If we have an element, but it's not a local element, then we
369  // either need to have a serialized vector or we need to find a
370  // local element sharing the same point.
371  if (element &&
372  (element->processor_id() != this->processor_id()) &&
373  _vector.type() != SERIAL)
374  {
375  // look for a local element containing the point
376  std::set<const Elem*> point_neighbors;
377  element->find_point_neighbors(p, point_neighbors);
378  element = NULL;
379  std::set<const Elem*>::const_iterator it = point_neighbors.begin();
380  const std::set<const Elem*>::const_iterator end = point_neighbors.end();
381  for (; it != end; ++it)
382  {
383  const Elem* elem = *it;
384  if (elem->processor_id() == this->processor_id())
385  {
386  element = elem;
387  break;
388  }
389  }
390  }
391 
392  if (!element)
393  {
394  output.resize(0);
395  }
396  else
397  {
398  // resize the output vector to the number of output values
399  // that the user told us
400  output.resize (this->_system_vars.size());
401 
402 
403  {
404  const unsigned int dim = this->_eqn_systems.get_mesh().mesh_dimension();
405 
406 
407  /*
408  * Get local coordinates to feed these into compute_data().
409  * Note that the fe_type can safely be used from the 0-variable,
410  * since the inverse mapping is the same for all FEFamilies
411  */
412  const Point mapped_point (FEInterface::inverse_map (dim,
413  this->_dof_map.variable_type(0),
414  element,
415  p));
416 
417  std::vector<Point> point_list (1, mapped_point);
418 
419  // loop over all vars
420  for (unsigned int index=0; index < this->_system_vars.size(); index++)
421  {
422  /*
423  * the data for this variable
424  */
425  const unsigned int var = _system_vars[index];
426  const FEType& fe_type = this->_dof_map.variable_type(var);
427 
428  AutoPtr<FEBase> point_fe (FEBase::build(dim, fe_type));
429  const std::vector<std::vector<RealGradient> >& dphi = point_fe->get_dphi();
430  point_fe->reinit(element, &point_list);
431 
432  // where the solution values for the var-th variable are stored
433  std::vector<dof_id_type> dof_indices;
434  this->_dof_map.dof_indices (element, dof_indices, var);
435 
436  // interpolate the solution
437  Gradient grad(0.);
438 
439  for (unsigned int i=0; i<dof_indices.size(); i++)
440  grad.add_scaled(dphi[i][0], this->_vector(dof_indices[i]));
441 
442  output[index] = grad;
443  }
444  }
445  }
446 
447  // all done
448  return;
449 }
450 
451 
452 
453 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
455  const Real,
456  std::vector<Tensor>& output)
457 {
458  libmesh_assert (this->initialized());
459 
460  /* Ensure that in the case of a master mesh function, the
461  out-of-mesh mode is enabled either for both or for none. This is
462  important because the out-of-mesh mode is also communicated to
463  the point locator. Since this is time consuming, enable it only
464  in debug mode. */
465 #ifdef DEBUG
466  if (this->_master != NULL)
467  {
468  const MeshFunction* master =
469  libmesh_cast_ptr<const MeshFunction*>(this->_master);
471  {
472  libMesh::err << "ERROR: If you use out-of-mesh-mode in connection with master mesh functions, you must enable out-of-mesh mode for both the master and the slave mesh function." << std::endl;
473  libmesh_error();
474  }
475  }
476 #endif
477 
478  // locate the point in the other mesh
479  const Elem* element = this->_point_locator->operator()(p);
480 
481  // If we have an element, but it's not a local element, then we
482  // either need to have a serialized vector or we need to find a
483  // local element sharing the same point.
484  if (element &&
485  (element->processor_id() != this->processor_id()) &&
486  _vector.type() != SERIAL)
487  {
488  // look for a local element containing the point
489  std::set<const Elem*> point_neighbors;
490  element->find_point_neighbors(p, point_neighbors);
491  element = NULL;
492  std::set<const Elem*>::const_iterator it = point_neighbors.begin();
493  const std::set<const Elem*>::const_iterator end = point_neighbors.end();
494  for (; it != end; ++it)
495  {
496  const Elem* elem = *it;
497  if (elem->processor_id() == this->processor_id())
498  {
499  element = elem;
500  break;
501  }
502  }
503  }
504 
505  if (!element)
506  {
507  output.resize(0);
508  }
509  else
510  {
511  // resize the output vector to the number of output values
512  // that the user told us
513  output.resize (this->_system_vars.size());
514 
515 
516  {
517  const unsigned int dim = this->_eqn_systems.get_mesh().mesh_dimension();
518 
519 
520  /*
521  * Get local coordinates to feed these into compute_data().
522  * Note that the fe_type can safely be used from the 0-variable,
523  * since the inverse mapping is the same for all FEFamilies
524  */
525  const Point mapped_point (FEInterface::inverse_map (dim,
526  this->_dof_map.variable_type(0),
527  element,
528  p));
529 
530  std::vector<Point> point_list (1, mapped_point);
531 
532  // loop over all vars
533  for (unsigned int index=0; index < this->_system_vars.size(); index++)
534  {
535  /*
536  * the data for this variable
537  */
538  const unsigned int var = _system_vars[index];
539  const FEType& fe_type = this->_dof_map.variable_type(var);
540 
541  AutoPtr<FEBase> point_fe (FEBase::build(dim, fe_type));
542  const std::vector<std::vector<RealTensor> >& d2phi =
543  point_fe->get_d2phi();
544  point_fe->reinit(element, &point_list);
545 
546  // where the solution values for the var-th variable are stored
547  std::vector<dof_id_type> dof_indices;
548  this->_dof_map.dof_indices (element, dof_indices, var);
549 
550  // interpolate the solution
551  Tensor hess;
552 
553  for (unsigned int i=0; i<dof_indices.size(); i++)
554  hess.add_scaled(d2phi[i][0], this->_vector(dof_indices[i]));
555 
556  output[index] = hess;
557  }
558  }
559  }
560 
561  // all done
562  return;
563 }
564 #endif
565 
566 
567 
569 {
570  libmesh_assert (this->initialized());
571  return *_point_locator;
572 }
573 
575 {
576  libmesh_assert (this->initialized());
578  _out_of_mesh_mode = true;
579  _out_of_mesh_value = value;
580 }
581 
583 {
584  DenseVector<Number> v(1);
585  v(0) = value;
586  this->enable_out_of_mesh_mode(v);
587 }
588 
590 {
591  libmesh_assert (this->initialized());
593  _out_of_mesh_mode = false;
594 }
595 
596 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:06 UTC

Hosted By:
SourceForge.net Logo