libMesh
trilinos_epetra_vector.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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 // C++ includes
19 #include <limits>
20 
21 // Local Includes
22 #include "libmesh/trilinos_epetra_vector.h"
23 
24 #ifdef LIBMESH_TRILINOS_HAVE_EPETRA
25 
26 #include "libmesh/dense_subvector.h"
27 #include "libmesh/dense_vector.h"
28 #include "libmesh/parallel.h"
29 #include "libmesh/trilinos_epetra_matrix.h"
30 #include "libmesh/utility.h"
31 
32 // Trilinos Includes
33 #include "libmesh/ignore_warnings.h"
34 #include <Epetra_LocalMap.h>
35 #include <Epetra_Comm.h>
36 #include <Epetra_Map.h>
37 #include <Epetra_BlockMap.h>
38 #include <Epetra_Import.h>
39 #include <Epetra_Export.h>
40 #include <Epetra_Util.h>
41 #include <Epetra_IntSerialDenseVector.h>
42 #include <Epetra_SerialDenseVector.h>
43 #include <Epetra_Vector.h>
44 #include "libmesh/restore_warnings.h"
45 
46 namespace libMesh
47 {
48 
49 template <typename T>
50 T EpetraVector<T>::sum () const
51 {
52  libmesh_assert(this->closed());
53 
54  const unsigned int nl = _vec->MyLength();
55 
56  T sum=0.0;
57 
58  T * values = _vec->Values();
59 
60  for (unsigned int i=0; i<nl; i++)
61  sum += values[i];
62 
63  this->comm().sum(sum);
64 
65  return sum;
66 }
67 
68 template <typename T>
69 Real EpetraVector<T>::l1_norm () const
70 {
71  libmesh_assert(this->closed());
72 
73  Real value;
74 
75  _vec->Norm1(&value);
76 
77  return value;
78 }
79 
80 template <typename T>
81 Real EpetraVector<T>::l2_norm () const
82 {
83  libmesh_assert(this->closed());
84 
85  Real value;
86 
87  _vec->Norm2(&value);
88 
89  return value;
90 }
91 
92 template <typename T>
93 Real EpetraVector<T>::linfty_norm () const
94 {
95  libmesh_assert(this->closed());
96 
97  Real value;
98 
99  _vec->NormInf(&value);
100 
101  return value;
102 }
103 
104 template <typename T>
105 NumericVector<T> &
106 EpetraVector<T>::operator += (const NumericVector<T> & v)
107 {
108  libmesh_assert(this->closed());
109 
110  this->add(1., v);
111 
112  return *this;
113 }
114 
115 
116 
117 template <typename T>
118 NumericVector<T> &
119 EpetraVector<T>::operator -= (const NumericVector<T> & v)
120 {
121  libmesh_assert(this->closed());
122 
123  this->add(-1., v);
124 
125  return *this;
126 }
127 
128 
129 template <typename T>
130 NumericVector<T> &
131 EpetraVector<T>::operator /= (NumericVector<T> & v)
132 {
133  libmesh_assert(this->closed());
134  libmesh_assert_equal_to(size(), v.size());
135 
136  EpetraVector<T> & v_vec = cast_ref<EpetraVector<T> &>(v);
137 
138  _vec->ReciprocalMultiply(1.0, *v_vec._vec, *_vec, 0.0);
139 
140  return *this;
141 }
142 
143 
144 
145 
146 template <typename T>
147 void EpetraVector<T>::set (const numeric_index_type i_in, const T value_in)
148 {
149  int i = static_cast<int> (i_in);
150  T value = value_in;
151 
152  libmesh_assert_less (i_in, this->size());
153 
154  ReplaceGlobalValues(1, &i, &value);
155 
156  this->_is_closed = false;
157 }
158 
159 
160 
161 template <typename T>
162 void EpetraVector<T>::reciprocal()
163 {
164  // The Epetra::reciprocal() function takes a constant reference to *another* vector,
165  // and fills _vec with its reciprocal. Does that mean we can't pass *_vec as the
166  // argument?
167  // _vec->reciprocal( *_vec );
168 
169  // Alternatively, compute the reciprocal by hand... see also the add(T) member that does this...
170  const unsigned int nl = _vec->MyLength();
171 
172  T * values = _vec->Values();
173 
174  for (unsigned int i=0; i<nl; i++)
175  {
176  // Don't divide by zero (maybe only check this in debug mode?)
177  if (std::abs(values[i]) < std::numeric_limits<T>::min())
178  libmesh_error_msg("Error, divide by zero in DistributedVector<T>::reciprocal()!");
179 
180  values[i] = 1. / values[i];
181  }
182 
183  // Leave the vector in a closed state...
184  this->close();
185 }
186 
187 
188 
189 template <typename T>
190 void EpetraVector<T>::conjugate()
191 {
192  // EPetra is real, rendering this a no-op.
193 }
194 
195 
196 
197 template <typename T>
198 void EpetraVector<T>::add (const numeric_index_type i_in, const T value_in)
199 {
200  int i = static_cast<int> (i_in);
201  T value = value_in;
202 
203  libmesh_assert_less (i_in, this->size());
204 
205  SumIntoGlobalValues(1, &i, &value);
206 
207  this->_is_closed = false;
208 }
209 
210 
211 
212 template <typename T>
213 void EpetraVector<T>::add_vector (const T * v,
214  const std::vector<numeric_index_type> & dof_indices)
215 {
216  libmesh_assert_equal_to (sizeof(numeric_index_type), sizeof(int));
217 
218  SumIntoGlobalValues (dof_indices.size(),
219  (int *) &dof_indices[0],
220  const_cast<T *>(v));
221 }
222 
223 
224 
225 // TODO: fill this in after creating an EpetraMatrix
226 template <typename T>
227 void EpetraVector<T>::add_vector (const NumericVector<T> & v_in,
228  const SparseMatrix<T> & A_in)
229 {
230  const EpetraVector<T> * v = cast_ptr<const EpetraVector<T> *>(&v_in);
231  const EpetraMatrix<T> * A = cast_ptr<const EpetraMatrix<T> *>(&A_in);
232 
233  // FIXME - does Trilinos let us do this *without* memory allocation?
234  UniquePtr<NumericVector<T>> temp = v->zero_clone();
235  EpetraVector<T> * temp_v = cast_ptr<EpetraVector<T> *>(temp.get());
236  A->mat()->Multiply(false, *v->_vec, *temp_v->_vec);
237  *this += *temp;
238 }
239 
240 
241 
242 // TODO: fill this in after creating an EpetraMatrix
243 template <typename T>
244 void EpetraVector<T>::add_vector_transpose (const NumericVector<T> & /* v_in */,
245  const SparseMatrix<T> & /* A_in */)
246 {
247  libmesh_not_implemented();
248 }
249 
250 
251 
252 template <typename T>
253 void EpetraVector<T>::add (const T v_in)
254 {
255  const unsigned int nl = _vec->MyLength();
256 
257  T * values = _vec->Values();
258 
259  for (unsigned int i=0; i<nl; i++)
260  values[i]+=v_in;
261 
262  this->_is_closed = false;
263 }
264 
265 
266 template <typename T>
267 void EpetraVector<T>::add (const NumericVector<T> & v)
268 {
269  this->add (1., v);
270 }
271 
272 
273 template <typename T>
274 void EpetraVector<T>::add (const T a_in, const NumericVector<T> & v_in)
275 {
276  const EpetraVector<T> * v = cast_ptr<const EpetraVector<T> *>(&v_in);
277 
278  libmesh_assert_equal_to (this->size(), v->size());
279 
280  _vec->Update(a_in,*v->_vec, 1.);
281 }
282 
283 
284 
285 template <typename T>
286 void EpetraVector<T>::insert (const T * v,
287  const std::vector<numeric_index_type> & dof_indices)
288 {
289  libmesh_assert_equal_to (sizeof(numeric_index_type), sizeof(int));
290 
291  ReplaceGlobalValues (dof_indices.size(),
292  (int *) &dof_indices[0],
293  const_cast<T *>(v));
294 }
295 
296 
297 
298 template <typename T>
299 void EpetraVector<T>::scale (const T factor_in)
300 {
301  _vec->Scale(factor_in);
302 }
303 
304 template <typename T>
306 {
307  _vec->Abs(*_vec);
308 }
309 
310 
311 template <typename T>
312 T EpetraVector<T>::dot (const NumericVector<T> & v_in) const
313 {
314  const EpetraVector<T> * v = cast_ptr<const EpetraVector<T> *>(&v_in);
315 
316  T result=0.0;
317 
318  _vec->Dot(*v->_vec, &result);
319 
320  return result;
321 }
322 
323 
324 template <typename T>
325 void EpetraVector<T>::pointwise_mult (const NumericVector<T> & vec1,
326  const NumericVector<T> & vec2)
327 {
328  const EpetraVector<T> * v1 = cast_ptr<const EpetraVector<T> *>(&vec1);
329  const EpetraVector<T> * v2 = cast_ptr<const EpetraVector<T> *>(&vec2);
330 
331  _vec->Multiply(1.0, *v1->_vec, *v2->_vec, 0.0);
332 }
333 
334 
335 template <typename T>
336 NumericVector<T> &
337 EpetraVector<T>::operator = (const T s_in)
338 {
339  _vec->PutScalar(s_in);
340 
341  return *this;
342 }
343 
344 
345 
346 template <typename T>
347 NumericVector<T> &
348 EpetraVector<T>::operator = (const NumericVector<T> & v_in)
349 {
350  const EpetraVector<T> * v = cast_ptr<const EpetraVector<T> *>(&v_in);
351 
352  *this = *v;
353 
354  return *this;
355 }
356 
357 
358 
359 template <typename T>
360 EpetraVector<T> &
361 EpetraVector<T>::operator = (const EpetraVector<T> & v)
362 {
363  (*_vec) = *v._vec;
364 
365  // FIXME - what about our communications data?
366 
367  return *this;
368 }
369 
370 
371 
372 template <typename T>
373 NumericVector<T> &
374 EpetraVector<T>::operator = (const std::vector<T> & v)
375 {
376  T * values = _vec->Values();
377 
382  if (this->size() == v.size())
383  {
384  const unsigned int nl=this->local_size();
385  const unsigned int fli=this->first_local_index();
386 
387  for (unsigned int i=0;i<nl;i++)
388  values[i]=v[fli+i];
389  }
390 
395  else
396  {
397  libmesh_assert_equal_to (v.size(), this->local_size());
398 
399  const unsigned int nl=this->local_size();
400 
401  for (unsigned int i=0;i<nl;i++)
402  values[i]=v[i];
403  }
404 
405  return *this;
406 }
407 
408 
409 
410 template <typename T>
411 void EpetraVector<T>::localize (NumericVector<T> & v_local_in) const
412 {
413  EpetraVector<T> * v_local = cast_ptr<EpetraVector<T> *>(&v_local_in);
414 
415  Epetra_Map rootMap = Epetra_Util::Create_Root_Map( *_map, -1);
416  v_local->_vec->ReplaceMap(rootMap);
417 
418  Epetra_Import importer(v_local->_vec->Map(), *_map);
419  v_local->_vec->Import(*_vec, importer, Insert);
420 }
421 
422 
423 
424 template <typename T>
425 void EpetraVector<T>::localize (NumericVector<T> & v_local_in,
426  const std::vector<numeric_index_type> & /* send_list */) const
427 {
428  // TODO: optimize to sync only the send list values
429  this->localize(v_local_in);
430 
431  // EpetraVector<T> * v_local =
432  // cast_ptr<EpetraVector<T> *>(&v_local_in);
433 
434  // libmesh_assert(this->_map.get());
435  // libmesh_assert(v_local->_map.get());
436  // libmesh_assert_equal_to (v_local->local_size(), this->size());
437  // libmesh_assert_less_equal (send_list.size(), v_local->size());
438 
439  // Epetra_Import importer (*v_local->_map, *this->_map);
440 
441  // v_local->_vec->Import (*this->_vec, importer, Insert);
442 }
443 
444 
445 
446 template <typename T>
447 void EpetraVector<T>::localize (std::vector<T> & v_local,
448  const std::vector<numeric_index_type> & indices) const
449 {
450  // Create a "replicated" map for importing values. This is
451  // equivalent to creating a general Epetra_Map with
452  // NumGlobalElements == NumMyElements.
453  Epetra_LocalMap import_map(static_cast<int>(indices.size()),
454  /*IndexBase=*/0,
455  _map->Comm());
456 
457  // Get a pointer to the list of global elements for the map, and set
458  // all the values from indices.
459  int * import_map_global_elements = import_map.MyGlobalElements();
460  for (int i=0; i<import_map.NumMyElements(); ++i)
461  import_map_global_elements[i] = indices[i];
462 
463  // Create a new EpetraVector to import values into.
464  Epetra_Vector import_vector(import_map);
465 
466  // Set up an "Import" object which associates the two maps.
467  Epetra_Import import_object(import_map, *_map);
468 
469  // Import the values
470  import_vector.Import(*_vec, import_object, Insert);
471 
472  // Get a pointer to the imported values array and the length of the
473  // array.
474  T * values = import_vector.Values();
475  int import_vector_length = import_vector.MyLength();
476 
477  // Copy the imported values into v_local
478  v_local.resize(import_vector_length);
479  for (int i=0; i<import_vector_length; ++i)
480  v_local[i] = values[i];
481 }
482 
483 
484 
485 template <typename T>
486 void EpetraVector<T>::localize (const numeric_index_type first_local_idx,
487  const numeric_index_type last_local_idx,
488  const std::vector<numeric_index_type> & send_list)
489 {
490  // Only good for serial vectors.
491  libmesh_assert_equal_to (this->size(), this->local_size());
492  libmesh_assert_greater (last_local_idx, first_local_idx);
493  libmesh_assert_less_equal (send_list.size(), this->size());
494  libmesh_assert_less (last_local_idx, this->size());
495 
496  const unsigned int my_size = this->size();
497  const unsigned int my_local_size = (last_local_idx - first_local_idx + 1);
498 
499  // Don't bother for serial cases
500  if ((first_local_idx == 0) &&
501  (my_local_size == my_size))
502  return;
503 
504  // Build a parallel vector, initialize it with the local
505  // parts of (*this)
506  EpetraVector<T> parallel_vec(this->comm(), PARALLEL);
507 
508  parallel_vec.init (my_size, my_local_size, true, PARALLEL);
509 
510  // Copy part of *this into the parallel_vec
511  for (numeric_index_type i=first_local_idx; i<=last_local_idx; i++)
512  parallel_vec.set(i,this->el(i));
513 
514  // localize like normal
515  parallel_vec.close();
516  parallel_vec.localize (*this, send_list);
517  this->close();
518 }
519 
520 
521 
522 template <typename T>
523 void EpetraVector<T>::localize (std::vector<T> & v_local) const
524 {
525  // This function must be run on all processors at once
526  parallel_object_only();
527 
528  const unsigned int n = this->size();
529  const unsigned int nl = this->local_size();
530 
531  libmesh_assert(this->_vec);
532 
533  v_local.clear();
534  v_local.reserve(n);
535 
536  // build up my local part
537  for (unsigned int i=0; i<nl; i++)
538  v_local.push_back((*this->_vec)[i]);
539 
540  this->comm().allgather (v_local);
541 }
542 
543 
544 
545 template <typename T>
546 void EpetraVector<T>::localize_to_one (std::vector<T> & v_local,
547  const processor_id_type pid) const
548 {
549  // This function must be run on all processors at once
550  parallel_object_only();
551 
552  const unsigned int n = this->size();
553  const unsigned int nl = this->local_size();
554 
555  libmesh_assert_less (pid, this->n_processors());
556  libmesh_assert(this->_vec);
557 
558  v_local.clear();
559  v_local.reserve(n);
560 
561 
562  // build up my local part
563  for (unsigned int i=0; i<nl; i++)
564  v_local.push_back((*this->_vec)[i]);
565 
566  this->comm().gather (pid, v_local);
567 }
568 
569 
570 
571 template <typename T>
572 void EpetraVector<T>::create_subvector(NumericVector<T> & /* subvector */,
573  const std::vector<numeric_index_type> & /* rows */) const
574 {
575  libmesh_not_implemented();
576 }
577 
578 
579 /*********************************************************************
580  * The following were copied (and slightly modified) from
581  * Epetra_FEVector.h in order to allow us to use a standard
582  * Epetra_Vector... which is more compatible with other Trilinos
583  * packages such as NOX. All of this code is originally under LGPL
584  *********************************************************************/
585 
586 //----------------------------------------------------------------------------
587 template <typename T>
588 int EpetraVector<T>::SumIntoGlobalValues(int numIDs,
589  const int * GIDs,
590  const double * values)
591 {
592  return( inputValues( numIDs, GIDs, values, true) );
593 }
594 
595 //----------------------------------------------------------------------------
596 template <typename T>
597 int EpetraVector<T>::SumIntoGlobalValues(const Epetra_IntSerialDenseVector & GIDs,
598  const Epetra_SerialDenseVector & values)
599 {
600  if (GIDs.Length() != values.Length()) {
601  return(-1);
602  }
603 
604  return( inputValues( GIDs.Length(), GIDs.Values(), values.Values(), true) );
605 }
606 
607 //----------------------------------------------------------------------------
608 template <typename T>
609 int EpetraVector<T>::SumIntoGlobalValues(int numIDs,
610  const int * GIDs,
611  const int * numValuesPerID,
612  const double * values)
613 {
614  return( inputValues( numIDs, GIDs, numValuesPerID, values, true) );
615 }
616 
617 //----------------------------------------------------------------------------
618 template <typename T>
619 int EpetraVector<T>::ReplaceGlobalValues(int numIDs,
620  const int * GIDs,
621  const double * values)
622 {
623  return( inputValues( numIDs, GIDs, values, false) );
624 }
625 
626 //----------------------------------------------------------------------------
627 template <typename T>
628 int EpetraVector<T>::ReplaceGlobalValues(const Epetra_IntSerialDenseVector & GIDs,
629  const Epetra_SerialDenseVector & values)
630 {
631  if (GIDs.Length() != values.Length()) {
632  return(-1);
633  }
634 
635  return( inputValues( GIDs.Length(), GIDs.Values(), values.Values(), false) );
636 }
637 
638 //----------------------------------------------------------------------------
639 template <typename T>
640 int EpetraVector<T>::ReplaceGlobalValues(int numIDs,
641  const int * GIDs,
642  const int * numValuesPerID,
643  const double * values)
644 {
645  return( inputValues( numIDs, GIDs, numValuesPerID, values, false) );
646 }
647 
648 //----------------------------------------------------------------------------
649 template <typename T>
650 int EpetraVector<T>::inputValues(int numIDs,
651  const int * GIDs,
652  const double * values,
653  bool accumulate)
654 {
655  if (accumulate) {
656  libmesh_assert(last_edit == 0 || last_edit == 2);
657  last_edit = 2;
658  } else {
659  libmesh_assert(last_edit == 0 || last_edit == 1);
660  last_edit = 1;
661  }
662 
663  //Important note!! This method assumes that there is only 1 point
664  //associated with each element.
665 
666  for (int i=0; i<numIDs; ++i) {
667  if (_vec->Map().MyGID(GIDs[i])) {
668  if (accumulate) {
669  _vec->SumIntoGlobalValue(GIDs[i], 0, 0, values[i]);
670  }
671  else {
672  _vec->ReplaceGlobalValue(GIDs[i], 0, 0, values[i]);
673  }
674  }
675  else {
676  if (!ignoreNonLocalEntries_) {
677  EPETRA_CHK_ERR( inputNonlocalValue(GIDs[i], values[i], accumulate) );
678  }
679  }
680  }
681 
682  return(0);
683 }
684 
685 //----------------------------------------------------------------------------
686 template <typename T>
687 int EpetraVector<T>::inputValues(int numIDs,
688  const int * GIDs,
689  const int * numValuesPerID,
690  const double * values,
691  bool accumulate)
692 {
693  if (accumulate) {
694  libmesh_assert(last_edit == 0 || last_edit == 2);
695  last_edit = 2;
696  } else {
697  libmesh_assert(last_edit == 0 || last_edit == 1);
698  last_edit = 1;
699  }
700 
701  int offset=0;
702  for (int i=0; i<numIDs; ++i) {
703  int numValues = numValuesPerID[i];
704  if (_vec->Map().MyGID(GIDs[i])) {
705  if (accumulate) {
706  for (int j=0; j<numValues; ++j) {
707  _vec->SumIntoGlobalValue(GIDs[i], j, 0, values[offset+j]);
708  }
709  }
710  else {
711  for (int j=0; j<numValues; ++j) {
712  _vec->ReplaceGlobalValue(GIDs[i], j, 0, values[offset+j]);
713  }
714  }
715  }
716  else {
717  if (!ignoreNonLocalEntries_) {
718  EPETRA_CHK_ERR( inputNonlocalValues(GIDs[i], numValues,
719  &(values[offset]), accumulate) );
720  }
721  }
722  offset += numValues;
723  }
724 
725  return(0);
726 }
727 
728 //----------------------------------------------------------------------------
729 template <typename T>
730 int EpetraVector<T>::inputNonlocalValue(int GID, double value, bool accumulate)
731 {
732  int insertPoint = -1;
733 
734  //find offset of GID in nonlocalIDs_
735  int offset = Epetra_Util_binary_search(GID, nonlocalIDs_, numNonlocalIDs_,
736  insertPoint);
737  if (offset >= 0) {
738  //if offset >= 0
739  // put value in nonlocalCoefs_[offset][0]
740 
741  if (accumulate) {
742  nonlocalCoefs_[offset][0] += value;
743  }
744  else {
745  nonlocalCoefs_[offset][0] = value;
746  }
747  }
748  else {
749  //else
750  // insert GID in nonlocalIDs_
751  // insert 1 in nonlocalElementSize_
752  // insert value in nonlocalCoefs_
753 
754  int tmp1 = numNonlocalIDs_;
755  int tmp2 = allocatedNonlocalLength_;
756  int tmp3 = allocatedNonlocalLength_;
757  EPETRA_CHK_ERR( Epetra_Util_insert(GID, insertPoint, nonlocalIDs_,
758  tmp1, tmp2) );
759  --tmp1;
760  EPETRA_CHK_ERR( Epetra_Util_insert(1, insertPoint, nonlocalElementSize_,
761  tmp1, tmp3) );
762  double * values = new double[1];
763  values[0] = value;
764  EPETRA_CHK_ERR( Epetra_Util_insert(values, insertPoint, nonlocalCoefs_,
765  numNonlocalIDs_, allocatedNonlocalLength_) );
766  }
767 
768  return(0);
769 }
770 
771 //----------------------------------------------------------------------------
772 template <typename T>
773 int EpetraVector<T>::inputNonlocalValues(int GID,
774  int numValues,
775  const double * values,
776  bool accumulate)
777 {
778  int insertPoint = -1;
779 
780  //find offset of GID in nonlocalIDs_
781  int offset = Epetra_Util_binary_search(GID, nonlocalIDs_, numNonlocalIDs_,
782  insertPoint);
783  if (offset >= 0) {
784  //if offset >= 0
785  // put value in nonlocalCoefs_[offset][0]
786 
787  if (numValues != nonlocalElementSize_[offset]) {
788  libMesh::err << "Epetra_FEVector ERROR: block-size for GID " << GID << " is "
789  << numValues<<" which doesn't match previously set block-size of "
790  << nonlocalElementSize_[offset] << std::endl;
791  return(-1);
792  }
793 
794  if (accumulate) {
795  for (int j=0; j<numValues; ++j) {
796  nonlocalCoefs_[offset][j] += values[j];
797  }
798  }
799  else {
800  for (int j=0; j<numValues; ++j) {
801  nonlocalCoefs_[offset][j] = values[j];
802  }
803  }
804  }
805  else {
806  //else
807  // insert GID in nonlocalIDs_
808  // insert numValues in nonlocalElementSize_
809  // insert values in nonlocalCoefs_
810 
811  int tmp1 = numNonlocalIDs_;
812  int tmp2 = allocatedNonlocalLength_;
813  int tmp3 = allocatedNonlocalLength_;
814  EPETRA_CHK_ERR( Epetra_Util_insert(GID, insertPoint, nonlocalIDs_,
815  tmp1, tmp2) );
816  --tmp1;
817  EPETRA_CHK_ERR( Epetra_Util_insert(numValues, insertPoint, nonlocalElementSize_,
818  tmp1, tmp3) );
819  double * newvalues = new double[numValues];
820  for (int j=0; j<numValues; ++j) {
821  newvalues[j] = values[j];
822  }
823  EPETRA_CHK_ERR( Epetra_Util_insert(newvalues, insertPoint, nonlocalCoefs_,
824  numNonlocalIDs_, allocatedNonlocalLength_) );
825  }
826 
827  return(0);
828 }
829 
830 //----------------------------------------------------------------------------
831 template <typename T>
832 int EpetraVector<T>::GlobalAssemble(Epetra_CombineMode mode)
833 {
834  //In this method we need to gather all the non-local (overlapping) data
835  //that's been input on each processor, into the (probably) non-overlapping
836  //distribution defined by the map that 'this' vector was constructed with.
837 
838  //We don't need to do anything if there's only one processor or if
839  //ignoreNonLocalEntries_ is true.
840  if (_vec->Map().Comm().NumProc() < 2 || ignoreNonLocalEntries_) {
841  return(0);
842  }
843 
844 
845 
846  //First build a map that describes the data in nonlocalIDs_/nonlocalCoefs_.
847  //We'll use the arbitrary distribution constructor of Map.
848 
849  Epetra_BlockMap sourceMap(-1, numNonlocalIDs_,
850  nonlocalIDs_, nonlocalElementSize_,
851  _vec->Map().IndexBase(), _vec->Map().Comm());
852 
853  //Now build a vector to hold our nonlocalCoefs_, and to act as the source-
854  //vector for our import operation.
855  Epetra_MultiVector nonlocalVector(sourceMap, 1);
856 
857  int i,j;
858  for (i=0; i<numNonlocalIDs_; ++i) {
859  for (j=0; j<nonlocalElementSize_[i]; ++j) {
860  nonlocalVector.ReplaceGlobalValue(nonlocalIDs_[i], j, 0,
861  nonlocalCoefs_[i][j]);
862  }
863  }
864 
865  Epetra_Export exporter(sourceMap, _vec->Map());
866 
867  EPETRA_CHK_ERR( _vec->Export(nonlocalVector, exporter, mode) );
868 
869  destroyNonlocalData();
870 
871  return(0);
872 }
873 
874 
875 //----------------------------------------------------------------------------
876 template <typename T>
877 void EpetraVector<T>::FEoperatorequals(const EpetraVector & source)
878 {
879  (*_vec) = *(source._vec);
880 
881  destroyNonlocalData();
882 
883  if (source.allocatedNonlocalLength_ > 0) {
884  allocatedNonlocalLength_ = source.allocatedNonlocalLength_;
885  numNonlocalIDs_ = source.numNonlocalIDs_;
886  nonlocalIDs_ = new int[allocatedNonlocalLength_];
887  nonlocalElementSize_ = new int[allocatedNonlocalLength_];
888  nonlocalCoefs_ = new double *[allocatedNonlocalLength_];
889  for (int i=0; i<numNonlocalIDs_; ++i) {
890  int elemSize = source.nonlocalElementSize_[i];
891  nonlocalCoefs_[i] = new double[elemSize];
892  nonlocalIDs_[i] = source.nonlocalIDs_[i];
893  nonlocalElementSize_[i] = elemSize;
894  for (int j=0; j<elemSize; ++j) {
895  nonlocalCoefs_[i][j] = source.nonlocalCoefs_[i][j];
896  }
897  }
898  }
899 }
900 
901 
902 //----------------------------------------------------------------------------
903 template <typename T>
904 void EpetraVector<T>::destroyNonlocalData()
905 {
906  if (allocatedNonlocalLength_ > 0) {
907  delete [] nonlocalIDs_;
908  delete [] nonlocalElementSize_;
909  nonlocalIDs_ = libmesh_nullptr;
910  nonlocalElementSize_ = libmesh_nullptr;
911  for (int i=0; i<numNonlocalIDs_; ++i) {
912  delete [] nonlocalCoefs_[i];
913  }
914  delete [] nonlocalCoefs_;
915  nonlocalCoefs_ = libmesh_nullptr;
916  numNonlocalIDs_ = 0;
917  allocatedNonlocalLength_ = 0;
918  }
919  return;
920 }
921 
922 
923 //------------------------------------------------------------------
924 // Explicit instantiations
925 template class EpetraVector<Number>;
926 
927 } // namespace libMesh
928 
929 #endif // LIBMESH_TRILINOS_HAVE_EPETRA
OStreamProxy err
bool closed()
Checks that the library has been closed.
Definition: libmesh.C:281
double abs(double a)
void scale(MeshBase &mesh, const Real xs, const Real ys=0., const Real zs=0.)
Scales the mesh.
uint8_t processor_id_type
Definition: id_types.h:99
const class libmesh_nullptr_t libmesh_nullptr
The libMesh namespace provides an interface to certain functionality in the library.
void sum(T &r, const Communicator &comm=Communicator_World)
libmesh_assert(j)
dof_id_type numeric_index_type
Definition: id_types.h:92
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static PetscErrorCode Mat * A
static const bool value
Definition: xdr_io.C:108
Iterator & operator=(const Iterator &rhs)
Assignment operator.
long double min(long double a, double b)
processor_id_type n_processors()
Definition: libmesh_base.h:88