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

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

Hosted By:
SourceForge.net Logo