Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
FEVector/ExecuteTestProblems.cpp
Go to the documentation of this file.
1//@HEADER
2// ************************************************************************
3//
4// Epetra: Linear Algebra Services Package
5// Copyright 2011 Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ************************************************************************
40//@HEADER
41
42
43#include "Epetra_BLAS.h"
44#include "ExecuteTestProblems.h"
45#include "Epetra_Comm.h"
46#include "Epetra_Vector.h"
49#include <vector>
50
51int MultiVectorTests(const Epetra_BlockMap & Map, int NumVectors, bool verbose)
52{
53 (void)NumVectors;
54 const Epetra_Comm & Comm = Map.Comm();
55 int ierr = 0;
56
57 /* get number of processors and the name of this processor */
58
59 // int NumProc = Comm.getNumProc();
60 int MyPID = Comm.MyPID();
61
62 // Construct FEVector
63
64 if (verbose&&MyPID==0) cout << "constructing Epetra_FEVector" << endl;
65
66 Epetra_FEVector A(Map, 1);
67
68 //For an extreme test, we'll have each processor sum-in a 1.0 for All
69 //global ids.
70
71 int minGID = Map.MinAllGID();
72 int numGlobalIDs = Map.NumGlobalElements();
73
74 //For now we're going to have just one point associated with
75 //each GID (element).
76
77 int* ptIndices = new int[numGlobalIDs];
78 double* ptCoefs = new double[numGlobalIDs];
79
80 Epetra_IntSerialDenseVector epetra_indices(View, ptIndices, numGlobalIDs);
81 Epetra_SerialDenseVector epetra_coefs(View, ptCoefs, numGlobalIDs);
82
83 {for(int i=0; i<numGlobalIDs; ++i) {
84 ptIndices[i] = minGID+i;
85 ptCoefs[i] = 1.0;
86 }}
87
88 if (verbose&&MyPID==0) {
89 cout << "calling A.SumIntoGlobalValues with " << numGlobalIDs << " values"<<endl;
90 }
91 EPETRA_TEST_ERR( A.SumIntoGlobalValues(numGlobalIDs, ptIndices, ptCoefs), ierr);
92
93 if (verbose&&MyPID==0) {
94 cout << "calling A.SumIntoGlobalValues with " << numGlobalIDs << " values"<<endl;
95 }
96 EPETRA_TEST_ERR( A.SumIntoGlobalValues(epetra_indices, epetra_coefs), ierr);
97
98 if (verbose&&MyPID==0) {
99 cout << "calling A.GlobalAssemble()" << endl;
100 }
101
102 EPETRA_TEST_ERR( A.GlobalAssemble(), ierr );
103
104 if (verbose&&MyPID==0) {
105 cout << "after globalAssemble"<<endl;
106 }
107 if (verbose) {
108 A.Print(cout);
109 }
110
111 //now do a quick test of the copy constructor
112 Epetra_FEVector B(A);
113
114 double nrm2a, nrm2b;
115 A.Norm2(&nrm2a);
116 B.Norm2(&nrm2b);
117
118 if (nrm2a != nrm2b) {
119 cerr << "copy-constructor test failed, norm of copy doesn't equal"
120 << " norm of original."<<endl;
121 return(-1);
122 }
123
124 delete [] ptIndices;
125 delete [] ptCoefs;
126
127 return(ierr);
128}
129
130int fevec0(Epetra_Comm& Comm, bool verbose)
131{
132 int ierr = 0;
133 int NumGlobalRows = 4;
134 int indexBase = 0;
135 Epetra_Map Map(NumGlobalRows, indexBase, Comm);
136
137 int Numprocs = Comm.NumProc();
138 int MyPID = Comm.MyPID();
139
140 if (Numprocs != 2) return(0);
141
142
143 int NumCols = 3;
144 int* Indices = new int[NumCols];
145
146 double* Values = new double[NumCols];
147
148// Create vectors
149
150 Epetra_FEVector b(Map, 1);
151 Epetra_FEVector x0(Map, 1);
152
153// source terms
154 NumCols = 2;
155
156 if(MyPID==0) // indices corresponding to element 0 on processor 0
157 {
158 Indices[0] = 0;
159 Indices[1] = 3;
160
161 Values[0] = 1./2.;
162 Values[1] = 1./2.;
163
164 }
165 else
166 {
167 Indices[0] = 1;
168 Indices[1] = 2;
169
170 Values[0] = 0;
171 Values[1] = 0;
172 }
173
174 EPETRA_TEST_ERR( b.SumIntoGlobalValues(NumCols, Indices, Values),
175 ierr);
176
178
179 if (verbose&&MyPID==0) {
180 cout << "b:"<<endl;
181 }
182
183 if (verbose) {
184 b.Print(cout);
185 }
186
187 x0 = b;
188
189 if (verbose&&MyPID==0) {
190 cout << "x:"<<endl;
191 }
192
193 if (verbose) {
194 x0.Print(cout);
195 }
196
197 delete [] Values;
198 delete [] Indices;
199
200 return(0);
201}
202
203int fevec1(Epetra_Comm& Comm, bool verbose)
204{
205 int Numprocs = Comm.NumProc();
206
207 if (Numprocs != 2) return(0);
208 int MyPID = Comm.MyPID();
209
210 int ierr = 0;
211 int NumGlobalRows = 6;
212 const int NumVectors = 4;
213 int indexBase = 0;
214 Epetra_Map Map(NumGlobalRows, indexBase, Comm);
215
216 const int Num = 4;
217 int Indices[Num];
218
219 double Values[Num];
220
221// Create vectors
222
223 Epetra_FEVector b(Map, NumVectors);
224 Epetra_FEVector x0(Map, NumVectors);
225
226// source terms
227
228 if(MyPID==0) // indices corresponding to element 0 on processor 0
229 {
230 Indices[0] = 0;
231 Indices[1] = 1;
232 Indices[2] = 4;
233 Indices[3] = 5;
234
235 Values[0] = 1./2.;
236 Values[1] = 1./2.;
237 Values[2] = 1./2.;
238 Values[3] = 1./2.;
239
240 }
241 else
242 {
243 Indices[0] = 1;
244 Indices[1] = 2;
245 Indices[2] = 3;
246 Indices[3] = 4;
247
248 Values[0] = 0;
249 Values[1] = 0;
250 Values[2] = 0;
251 Values[3] = 0;
252 }
253
254 for(int i=0; i<NumVectors; ++i) {
255 EPETRA_TEST_ERR( b.SumIntoGlobalValues(Num, Indices, Values, i),
256 ierr);
257 }
258
260
261 double nrm2[NumVectors];
262
263 b.Norm2(nrm2);
264
265 for(int i=1; i<NumVectors; ++i) {
266 if (fabs(nrm2[i]-nrm2[0]) > 1.e-12) {
267 EPETRA_TEST_ERR(-1, ierr);
268 return(-1);
269 }
270 }
271
272
273 //now sum-in again, to make sure the previous call to GlobalAssemble
274 //didn't do something nasty to internal non-local data structures.
275 //(This is a specific case that has bitten me. Hence this test...)
276 for(int i=0; i<NumVectors; ++i) {
277 EPETRA_TEST_ERR( b.SumIntoGlobalValues(Num, Indices, Values, i),
278 ierr);
279 }
280
281 //and now GlobalAssemble again...
283
284
285 if (verbose&&MyPID==0) {
286 cout << "b:"<<endl;
287 }
288
289 if (verbose) {
290 b.Print(cout);
291 }
292
293 x0 = b;
294
295 if (verbose&&MyPID==0) {
296 cout << "x:"<<endl;
297 }
298
299 if (verbose) {
300 x0.Print(cout);
301 }
302
303 return(0);
304}
305
306int fevec2(Epetra_Comm& Comm, bool verbose)
307{
308 int ierr = 0;
309 int NumGlobalElems = 4;
310 int elemSize = 3;
311 int indexBase = 0;
312 Epetra_BlockMap Map(NumGlobalElems, elemSize, indexBase, Comm);
313
314 int Numprocs = Comm.NumProc();
315 int MyPID = Comm.MyPID();
316
317 if (Numprocs != 2) return(0);
318
319 int NumCols = 3;
320 int* Indices = new int[NumCols];
321 int* numValuesPerID = new int[NumCols];
322 for(int i=0; i<NumCols; ++i) {
323 numValuesPerID[i] = elemSize;
324 }
325
326 double* Values = new double[NumCols*elemSize];
327
328// Create vectors
329
330 Epetra_FEVector b(Map, 1);
331 Epetra_FEVector x0(Map, 1);
332
333// source terms
334 NumCols = 2;
335
336 if(MyPID==0) // indices corresponding to element 0 on processor 0
337 {
338 Indices[0] = 0;
339 Indices[1] = 3;
340
341 Values[0] = 1./2.;
342 Values[1] = 1./2.;
343 Values[2] = 1./2.;
344 Values[3] = 1./2.;
345 Values[4] = 1./2.;
346 Values[5] = 1./2.;
347
348 }
349 else
350 {
351 Indices[0] = 1;
352 Indices[1] = 2;
353
354 Values[0] = 0;
355 Values[1] = 0;
356 Values[2] = 0;
357 Values[3] = 0;
358 Values[4] = 0;
359 Values[5] = 0;
360 }
361
362 EPETRA_TEST_ERR( b.SumIntoGlobalValues(NumCols, Indices,
363 numValuesPerID, Values),
364 ierr);
365
367
368 if (verbose&&MyPID==0) {
369 cout << "b:"<<endl;
370 }
371
372 if (verbose) {
373 b.Print(cout);
374 }
375
376 x0 = b;
377
378 if (verbose&&MyPID==0) {
379 cout << "x:"<<endl;
380 }
381
382 if (verbose) {
383 x0.Print(cout);
384 }
385
386 delete [] Values;
387 delete [] Indices;
388 delete [] numValuesPerID;
389
390 return(0);
391}
392
393int fevec3(Epetra_Comm& Comm, bool verbose)
394{
395 int ierr = 0;
396 int NumGlobalElems = 4;
397 int elemSize = 40;
398 int indexBase = 0;
399 Epetra_BlockMap Map(NumGlobalElems, elemSize, indexBase, Comm);
400
401 int Numprocs = Comm.NumProc();
402 int MyPID = Comm.MyPID();
403
404 if (Numprocs != 2) return(0);
405
406 int NumCols = 3;
407 int* Indices = new int[NumCols];
408 int* numValuesPerID = new int[NumCols];
409 for(int i=0; i<NumCols; ++i) {
410 numValuesPerID[i] = elemSize;
411 }
412
413 double* Values = new double[NumCols*elemSize];
414
415// Create vectors
416
417 Epetra_FEVector b(Map, 1);
418 Epetra_FEVector x0(Map, 1);
419
420// source terms
421 NumCols = 2;
422
423 if(MyPID==0) // indices corresponding to element 0 on processor 0
424 {
425 Indices[0] = 0;
426 Indices[1] = 3;
427
428 for(int ii=0; ii<NumCols*elemSize; ++ii) {
429 Values[ii] = 1./2.;
430 }
431
432 }
433 else
434 {
435 Indices[0] = 1;
436 Indices[1] = 2;
437
438 for(int ii=0; ii<NumCols*elemSize; ++ii) {
439 Values[ii] = 0.;
440 }
441
442 }
443
444 EPETRA_TEST_ERR( b.SumIntoGlobalValues(NumCols, Indices,
445 numValuesPerID, Values),
446 ierr);
447
449
450 if (verbose&&MyPID==0) {
451 cout << "b:"<<endl;
452 }
453
454 if (verbose) {
455 b.Print(cout);
456 }
457
458 x0 = b;
459
460 if (verbose&&MyPID==0) {
461 cout << "x:"<<endl;
462 }
463
464 if (verbose) {
465 x0.Print(cout);
466 }
467
468 delete [] Values;
469 delete [] Indices;
470 delete [] numValuesPerID;
471
472 return(0);
473}
474
475int fevec4(Epetra_Comm& Comm, bool verbose)
476{
477 int NumElements = 4;
478 Epetra_Map Map(NumElements, 0, Comm);
479 Epetra_FEVector x1(Map);
480 const double value = 1.;
481 x1.PutScalar (value);
482 // replace one element by itself. processor 0
483 // does not own this element
484 const int GID = 3;
485 x1.ReplaceGlobalValues(1, &GID, &value);
487
488 if (Map.MyGID(3)) {
489 //insist that the value for GID==3 is 1:
490 if (std::abs(x1.Values()[Map.LID(3)] - 1) > 1.e-9) return -1;
491 }
492
493 std::cout << x1;
494
495 Comm.Barrier();
496
497 // re-apply GlobalAssemble. Nothing should
498 // happen
500 std::cout << x1;
501 if (Map.MyGID(3)) {
502 //insist that the value for GID==3 is 1:
503 if (std::abs(x1.Values()[Map.LID(3)] - 1) > 1.e-9) return -1;
504 }
505
506 return 0;
507}
508
509int fevec5(Epetra_Comm& Comm, bool verbose)
510{
511 int NumElements = 4;
512 Epetra_Map Map(NumElements, 0, Comm);
513 Epetra_FEVector x1(Map);
514 x1.PutScalar (0);
515
516 // let all processors set global entry 0 to 1
517 const int GID = 0;
518 const double value = 1;
519 x1.ReplaceGlobalValues(1, &GID, &value);
521 if (Comm.MyPID()==0)
522 std::cout << "Entry " << GID << " after construct & set: "
523 << x1[0][0] << std::endl;
524
525 // copy vector
526 Epetra_FEVector x2 (x1);
527
528 x2.PutScalar(0);
529
530 // re-apply 1 to the vector, but only on the
531 // owning processor. should be enough to set
532 // the value (as non-local data in x1 should
533 // have been eliminated after calling
534 // GlobalAssemble).
535 if (Comm.MyPID()==0)
536 x2.ReplaceGlobalValues(1, &GID, &value);
538
539 if (Comm.MyPID()==0)
540 std::cout << "Entry " << GID << " after copy & set: "
541 << x2[0][0] << std::endl;
542
543 return 0;
544}
545
546int fevec6(Epetra_Comm& Comm, bool verbose)
547{
548 int NumElements = 4;
549 Epetra_Map Map(NumElements, 0, Comm);
550 Epetra_FEVector x1(Map);
551 x1.PutScalar (0);
552
553 // let all processors set global entry 0 to 1
554 const int GID = 0;
555 const double value = 1;
556 x1.ReplaceGlobalValues(1, &GID, &value);
558 if (Comm.MyPID()==0)
559 std::cout << "Entry " << GID << " after construct & set: "
560 << x1[0][0] << std::endl;
561
562 x1.PutScalar(0);
563
564 // re-apply 1 to the vector, but only on the
565 // owning processor. should be enough to set
566 // the value (as non-local data in x1 should
567 // have been eliminated after calling
568 // GlobalAssemble).
569 if (Comm.MyPID()==0)
570 x1.ReplaceGlobalValues(1, &GID, &value);
572
573 if (Comm.MyPID()==0) {
574 std::cout << "Entry " << GID << " after PutScalar & set: "
575 << x1[0][0] << std::endl;
576 if (x1[0][0] != value) return -1;
577 }
578
579 return 0;
580}
581
582int fevec7(Epetra_Comm& Comm, bool verbose)
583{
584 const int NumVectors = 4;
585 const int NumElements = 4;
586 Epetra_Map Map(NumElements, 0, Comm);
587 std::vector<double> mydata(NumElements*NumVectors, 1.0);
588 Epetra_FEVector x1(View, Map, &mydata[0], NumElements, NumVectors);
589
590 x1.PutScalar (0);
591
592 // let all processors set global entry 0 to 1
593 const int GID = 0;
594 const double value = 1;
595 x1.ReplaceGlobalValues(1, &GID, &value);
597
598 if (Comm.MyPID()==0 && x1[0][0] != value) return -1;
599 return 0;
600}
601
int fevec3(Epetra_Comm &Comm, bool verbose)
int fevec0(Epetra_Comm &Comm, bool verbose)
int fevec2(Epetra_Comm &Comm, bool verbose)
int fevec6(Epetra_Comm &Comm, bool verbose)
int fevec4(Epetra_Comm &Comm, bool verbose)
int MultiVectorTests(const Epetra_BlockMap &Map, int NumVectors, bool verbose)
int fevec1(Epetra_Comm &Comm, bool verbose)
int fevec7(Epetra_Comm &Comm, bool verbose)
int fevec5(Epetra_Comm &Comm, bool verbose)
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
int MinAllGID() const
Returns the minimum global ID across the entire map.
int LID(int GID) const
Returns local ID of global ID, return -1 if not found on this processor.
int NumGlobalElements() const
Number of elements across all processors.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
bool MyGID(int GID_in) const
Returns true if the GID passed in belongs to the calling processor in this map, otherwise returns fal...
Epetra_Comm: The Epetra Communication Abstract Base Class.
Definition Epetra_Comm.h:73
virtual int NumProc() const =0
Returns total number of processes.
virtual int MyPID() const =0
Return my process ID.
virtual void Barrier() const =0
Epetra_Comm Barrier function.
Epetra Finite-Element Vector.
int SumIntoGlobalValues(int numIDs, const int *GIDs, const double *values, int vectorIndex=0)
Accumulate values into the vector, adding them to any values that already exist for the specified ind...
int GlobalAssemble(Epetra_CombineMode mode=Add, bool reuse_map_and_exporter=false)
Gather any overlapping/shared data into the non-overlapping partitioning defined by the Map that was ...
int ReplaceGlobalValues(int numIDs, const int *GIDs, const double *values, int vectorIndex=0)
Copy values into the vector overwriting any values that already exist for the specified indices.
Epetra_IntSerialDenseVector: A class for constructing and using dense vectors.
Epetra_Map: A class for partitioning vectors and matrices.
Definition Epetra_Map.h:119
virtual void Print(std::ostream &os) const
Print method.
double * Values() const
Get pointer to MultiVector values.
int Norm2(double *Result) const
Compute 2-norm of each vector in multi-vector.
int PutScalar(double ScalarConstant)
Initialize all values in a multi-vector with constant value.
Epetra_SerialDenseVector: A class for constructing and using dense vectors.
#define EPETRA_TEST_ERR(a, b)