EpetraExt Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
GLpApp_GLpYUEpetraDataPool.cpp
Go to the documentation of this file.
1/*
2//@HEADER
3// ***********************************************************************
4//
5// EpetraExt: Epetra Extended - Linear Algebra Services Package
6// Copyright (2011) Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39//
40// ***********************************************************************
41//@HEADER
42*/
43
44
45#include <cstdlib>
46#include <algorithm>
47#include <iostream>
48
51
52#include "Amesos_Klu.h"
56#include "Epetra_BLAS.h"
57#include "Epetra_CrsMatrix.h"
58#include "Epetra_Export.h"
59#include "Epetra_FECrsMatrix.h"
60#include "Epetra_FEVector.h"
61#include "Epetra_Import.h"
63#include "Epetra_LAPACK.h"
64#include "Epetra_Map.h"
65#include "Epetra_MultiVector.h"
67#include "Epetra_Vector.h"
68#include "Teuchos_ParameterList.hpp"
69#include "Teuchos_RCP.hpp"
70#include "Teuchos_Assert.hpp"
71#include "Teuchos_VerboseObject.hpp"
72
73#ifdef HAVE_MPI
74# include "Epetra_MpiComm.h"
75# include "mpi.h"
76#else
77# include "Epetra_SerialComm.h"
78#endif
79
80
81// 2008/09/04: Added to address failed tests (see bug 4040)
82using namespace std;
83
84
85// Define to see all debug output for mesh generation
86//#define GLPYUEPETRA_DATAPOOL_DUMP_ALL_MESH
87
88
89namespace GLpApp {
90
91//
92// Declarations
93//
94
95const double GLp_pi = 3.14159265358979323846;
96
97std::ostream& operator <<(std::ostream &, const Usr_Par &);
98
99bool CrsMatrix2MATLAB(const Epetra_CrsMatrix &, std::ostream &);
100
101bool Vector2MATLAB( const Epetra_Vector &, std::ostream &);
102
103bool FEVector2MATLAB( const Epetra_FEVector &, std::ostream &);
104
105int quadrature(const int, const int, Epetra_SerialDenseMatrix &,
107
108int meshreader(
109 const Epetra_Comm &,
116 const char geomFileBase[],
117 const bool trace = true,
118 const bool dumpAll = false
119 );
120
126 const Usr_Par &,
129
137 Teuchos::RCP<Epetra_FECrsMatrix> &,
138 Teuchos::RCP<Epetra_FEVector> &);
139
147 Teuchos::RCP<Epetra_FECrsMatrix> &,
148 Teuchos::RCP<Epetra_FEVector> &,
149 bool);
150
158 Teuchos::RCP<Epetra_FECrsMatrix> &,
159 Teuchos::RCP<Epetra_FECrsMatrix> &,
160 Teuchos::RCP<Epetra_FEVector> &);
161
162int assemble_bdry(
163 const Epetra_Comm &Comm
164 ,const Epetra_IntSerialDenseVector &ipindx
165 ,const Epetra_SerialDenseMatrix &ipcoords
166 ,const Epetra_IntSerialDenseVector &pindx
167 ,const Epetra_SerialDenseMatrix &pcoords
170 ,Teuchos::RCP<Epetra_FECrsMatrix> *B
171 ,Teuchos::RCP<Epetra_FECrsMatrix> *R
172 );
173
174int nonlinvec(const Epetra_Comm &,
180 const Teuchos::RCP<const Epetra_MultiVector> &,
181 Teuchos::RCP<Epetra_FEVector> &);
182
183int nonlinjac(const Epetra_Comm &,
189 const Teuchos::RCP<const Epetra_MultiVector> &,
190 Teuchos::RCP<Epetra_FECrsMatrix> &);
191
192int nonlinhessvec(const Epetra_Comm &,
198 const Teuchos::RCP<const Epetra_MultiVector> &,
199 const Teuchos::RCP<const Epetra_MultiVector> &,
200 const Teuchos::RCP<const Epetra_MultiVector> &,
201 Teuchos::RCP<Epetra_FEVector> &);
202
203int compproduct(Epetra_SerialDenseVector &, double *, double *);
204
205int compproduct(Epetra_SerialDenseVector &, double *, double *, double *);
206
208
210
211int quadrature(
212 const int, const int, Epetra_SerialDenseMatrix &,
214
216
218
220
221//
222// GLpYUEpetraDataPool
223//
224
226 Teuchos::RCP<const Epetra_Comm> const& commptr
227 ,const double beta
228 ,const double len_x
229 ,const double len_y
230 ,const int local_nx
231 ,const int local_ny
232 ,const char myfile[]
233 ,const bool trace
234 )
235 :commptr_(commptr)
236 ,beta_(beta)
237{
238 ipcoords_ = Teuchos::rcp( new Epetra_SerialDenseMatrix() );
239 ipindx_ = Teuchos::rcp( new Epetra_IntSerialDenseVector() );
240 pcoords_ = Teuchos::rcp( new Epetra_SerialDenseMatrix() );
241 pindx_ = Teuchos::rcp( new Epetra_IntSerialDenseVector() );
242 t_ = Teuchos::rcp( new Epetra_IntSerialDenseMatrix() );
243 e_ = Teuchos::rcp( new Epetra_IntSerialDenseMatrix() );
244
245 if( myfile && myfile[0] != '\0' ) {
247 }
248 else {
250 commptr_->NumProc(),commptr_->MyPID()
251 ,len_x,len_y,local_nx,local_ny,2 // 2==Neuman boundary conditions!
252 ,&*ipindx_,&*ipcoords_,&*pindx_,&*pcoords_,&*t_,&*e_
253#ifdef GLPYUEPETRA_DATAPOOL_DUMP_ALL_MESH
254 ,&*Teuchos::VerboseObjectBase::getDefaultOStream(),true
255#else
256 ,NULL,false
257#endif
258 );
259 }
260
261 // Assemble volume and boundary mass and stiffness matrices, and the right-hand side of the PDE.
262 assemble( *commptr, *ipindx_, *ipcoords_, *pindx_, *pcoords_, *t_, *e_, A_, H_, b_ );
263 assemble_bdry( *commptr, *ipindx_, *ipcoords_, *pindx_, *pcoords_, *t_, *e_, &B_, &R_ );
264
265 // Set desired state q.
266 Epetra_Map standardmap(A_->DomainMap());
267 q_ = Teuchos::rcp(new Epetra_FEVector(standardmap,1));
268 int * qintvalues = new int[standardmap.NumMyElements()];
269 double * qdvalues = new double[standardmap.NumMyElements()];
270 standardmap.MyGlobalElements(qintvalues);
271 for (int i = 0; i < standardmap.NumMyElements(); i++)
272 qdvalues[i]=cos( GLp_pi* ((*ipcoords_)(i,0)) ) * cos( GLp_pi* ((*ipcoords_)(i,1)) );
273 q_->ReplaceGlobalValues(standardmap.NumMyElements(), qintvalues, qdvalues);
274 q_->GlobalAssemble();
275}
276
278{
279
280 // Dynamic cast back to Epetra vectors here.
281 Teuchos::RCP<const Epetra_MultiVector> ey =
282 (Teuchos::dyn_cast<GenSQP::YUEpetraVector>(const_cast<GenSQP::Vector&>(x))).getYVector();
283
284 computeNy(ey);
285
286 computeNpy(ey);
287
289
290}
291
293 const Teuchos::RCP<const Epetra_MultiVector> & rhsy,
294 const Teuchos::RCP<const Epetra_MultiVector> & rhsu,
295 const Teuchos::RCP<const Epetra_MultiVector> & rhsp,
296 const Teuchos::RCP<Epetra_MultiVector> & y,
297 const Teuchos::RCP<Epetra_MultiVector> & u,
298 const Teuchos::RCP<Epetra_MultiVector> & p,
299 double tol )
300{
301/*
302 int systemChoice = 1; // 1 for full KKT system solve, 2 for Schur complement solve
303 int solverChoice = 12; // These options are for the full KKT system solve.
304 // 11 for AztecOO with built-in Schwarz DD preconditioning and ILU on subdomains
305 // 12 for AztecOO with IFPACK Schwarz DD preconditioning and Umfpack on subdomains
306 // 13 for a direct sparse solver (Umfpack, KLU)
307
308 if (systemChoice == 1) {
309 // We're using the full KKT system formulation to solve the augmented system.
310
311 Epetra_Map standardmap(A_->DomainMap());
312 int numstates = standardmap.NumGlobalElements();
313 Epetra_Map bdryctrlmap(B_->DomainMap());
314 int numcontrols = bdryctrlmap.NumGlobalElements();
315 Epetra_Vector rhs( (Epetra_BlockMap&)Augmat_->RangeMap() );
316 Epetra_Vector soln( (Epetra_BlockMap&)Augmat_->RangeMap() );
317 soln.Random();
318
319 std::vector<double> values(rhsy->MyLength() + rhsu->MyLength() + rhsp->MyLength());
320 std::vector<int> indices(rhsy->MyLength() + rhsu->MyLength() + rhsp->MyLength());
321 ((Epetra_BlockMap&)Augmat_->RangeMap()).MyGlobalElements(&indices[0]);
322
323 for (int i=0; i<rhsy->MyLength(); i++) {
324 values[i] = (*((*rhsy)(0)))[i];
325 }
326 for (int i=0; i<rhsu->MyLength(); i++) {
327 values[i+rhsy->MyLength()] = (*((*rhsu)(0)))[i];
328 }
329 for (int i=0; i<rhsp->MyLength(); i++) {
330 values[i+rhsy->MyLength()+rhsu->MyLength()] = (*((*rhsp)(0)))[i];
331 }
332
333 rhs.ReplaceGlobalValues(rhsy->MyLength() + rhsu->MyLength() + rhsp->MyLength(), &values[0], &indices[0]);
334
335 if (solverChoice == 11) {
336 int Overlap = 3;
337 int ival = 4;
338
339 AztecOO::AztecOO kktsolver(&(*Augmat_), &soln, &rhs);
340 kktsolver.SetAztecOption( AZ_solver, AZ_gmres );
341 kktsolver.SetAztecOption( AZ_precond, AZ_dom_decomp );
342 kktsolver.SetAztecOption(AZ_subdomain_solve, AZ_ilu);
343 //kktsolver.SetAztecOption( AZ_kspace, 2*numstates+numcontrols );
344 kktsolver.SetAztecOption( AZ_kspace, 9000 );
345 kktsolver.SetAztecOption(AZ_overlap,Overlap);
346 kktsolver.SetAztecOption(AZ_graph_fill,ival);
347 //kktsolver.SetAztecOption(AZ_poly_ord, ival);
348 //kktsolver.SetAztecParam(AZ_drop, 1e-9);
349 kktsolver.SetAztecParam(AZ_athresh, 1e-5);
350 //kktsolver.SetAztecParam(AZ_rthresh, 0.0);
351 kktsolver.SetAztecOption( AZ_reorder, 0 );
352 //kktsolver.SetAztecParam44( AZ_ilut_fill, 1.5 );
353 kktsolver.SetAztecOption( AZ_output, AZ_last );
354 //kktsolver.Iterate(2*numstates+numcontrols,1e-12);
355 kktsolver.Iterate(9000,1e-11);
356 //std::cout << soln;
357 }
358 else if (solverChoice == 12) {
359 // =============================================================== //
360 // B E G I N N I N G O F I F P A C K C O N S T R U C T I O N //
361 // =============================================================== //
362
363 Teuchos::ParameterList List;
364
365 // allocates an IFPACK factory. No data is associated
366 // to this object (only method Create()).
367 Ifpack Factory;
368
369 // create the preconditioner. For valid PrecType values,
370 // please check the documentation
371 std::string PrecType = "Amesos";
372 int OverlapLevel = 2; // must be >= 0. If Comm.NumProc() == 1,
373 // it is ignored.
374
375 Ifpack_Preconditioner* Prec = Factory.Create(PrecType, &(*Augmat_), OverlapLevel);
376 assert(Prec != 0);
377
378 // specify the Amesos solver to be used.
379 // If the selected solver is not available,
380 // IFPACK will try to use Amesos' KLU (which is usually always
381 // compiled). Amesos' serial solvers are:
382 // "Amesos_Klu", "Amesos_Umfpack", "Amesos_Superlu"
383 List.set("amesos: solver type", "Amesos_Umfpack");
384
385 // sets the parameters
386 IFPACK_CHK_ERR(Prec->SetParameters(List));
387
388 // initialize the preconditioner. At this point the matrix must
389 // have been FillComplete()'d, but actual values are ignored.
390 // At this call, Amesos will perform the symbolic factorization.
391 IFPACK_CHK_ERR(Prec->Initialize());
392
393 // Builds the preconditioners, by looking for the values of
394 // the matrix. At this call, Amesos will perform the
395 // numeric factorization.
396 IFPACK_CHK_ERR(Prec->Compute());
397
398 // =================================================== //
399 // E N D O F I F P A C K C O N S T R U C T I O N //
400 // =================================================== //
401
402 // need an Epetra_LinearProblem to define AztecOO solver
403 Epetra_LinearProblem Problem;
404 Problem.SetOperator(&(*Augmat_));
405 Problem.SetLHS(&soln);
406 Problem.SetRHS(&rhs);
407
408 // now we can allocate the AztecOO solver
409 AztecOO kktsolver(Problem);
410
411 // specify solver
412 kktsolver.SetAztecOption(AZ_solver,AZ_gmres);
413 kktsolver.SetAztecOption(AZ_kspace, 300 );
414 kktsolver.SetAztecOption(AZ_output,AZ_last);
415
416 // HERE WE SET THE IFPACK PRECONDITIONER
417 kktsolver.SetPrecOperator(Prec);
418
419 // .. and here we solve
420 kktsolver.Iterate(300,1e-12);
421
422 // delete the preconditioner
423 delete Prec;
424 }
425 else if (solverChoice == 13) {
426 Epetra_LinearProblem Problem;
427 Problem.SetOperator(&(*Augmat_));
428 Problem.SetLHS(&soln);
429 Problem.SetRHS(&rhs);
430
431 EpetraExt::LinearProblem_Reindex reindex(NULL);
432 Epetra_LinearProblem newProblem = reindex(Problem);
433
434 Amesos_Klu kktsolver(newProblem);
435
436 AMESOS_CHK_ERR(kktsolver.SymbolicFactorization());
437 AMESOS_CHK_ERR(kktsolver.NumericFactorization());
438 AMESOS_CHK_ERR(kktsolver.Solve());
439 kktsolver.PrintTiming();
440 }
441
442
443 for (int i=0; i<rhsy->MyLength(); i++) {
444 (*((*y)(0)))[i] = soln[i];
445 }
446 for (int i=0; i<rhsu->MyLength(); i++) {
447 (*((*u)(0)))[i] = soln[i+rhsy->MyLength()];
448 }
449 for (int i=0; i<rhsp->MyLength(); i++) {
450 (*((*p)(0)))[i] = soln[i+rhsy->MyLength()+rhsu->MyLength()];
451 }
452
453 }
454 else if (systemChoice == 2) {
455 // We're using the Schur complement formulation to solve the augmented system.
456
457 // Form linear operator.
458 GLpApp::SchurOp schurop(A_, B_, Npy_);
459
460 // Form Schur complement right-hand side.
461 Epetra_MultiVector ny( (Epetra_BlockMap&)Npy_->RangeMap(), 1);
462 Epetra_MultiVector ay( (Epetra_BlockMap&)A_->RangeMap(), 1);
463 Epetra_MultiVector schurrhs( (Epetra_BlockMap&)A_->RangeMap(), 1);
464 Epetra_MultiVector bu( (Epetra_BlockMap&)B_->RangeMap(), 1);
465 A_->Multiply(false, *rhsy, ay);
466 Npy_->Multiply(false, *rhsy, ny);
467 B_->Multiply(false, *rhsu, bu);
468 schurrhs.Update(1.0, ny, 1.0, ay, 0.0);
469 schurrhs.Update(-1.0, *rhsp, 1.0, bu, 1.0);
470
471 p->PutScalar(0.0);
472 Epetra_LinearProblem linprob(&schurop, &(*p), &schurrhs);
473 AztecOO::AztecOO Solver(linprob);
474 Solver.SetAztecOption( AZ_solver, AZ_cg );
475 Solver.SetAztecOption( AZ_precond, AZ_none );
476 Solver.SetAztecOption( AZ_output, AZ_none );
477 Solver.Iterate(8000,tol);
478
479 Epetra_MultiVector bp( (Epetra_BlockMap&)B_->DomainMap(), 1);
480 B_->Multiply(true, *p, bp);
481 u->Update(1.0, *rhsu, -1.0, bp, 0.0);
482
483 Epetra_MultiVector ap( (Epetra_BlockMap&)A_->DomainMap(), 1);
484 Epetra_MultiVector np( (Epetra_BlockMap&)A_->DomainMap(), 1);
485 A_->Multiply(true, *p, ap);
486 Npy_->Multiply(true, *p, np);
487 y->Update(1.0, *rhsy,0.0);
488 y->Update(-1.0, ap, -1.0, np, 1.0);
489 }
490*/
491 TEUCHOS_TEST_FOR_EXCEPT(true);
492 return 0;
493}
494
495Teuchos::RCP<const Epetra_Comm> GLpYUEpetraDataPool::getCommPtr() { return commptr_; }
496
497Teuchos::RCP<Epetra_FECrsMatrix> GLpYUEpetraDataPool::getA() { return A_; }
498
499Teuchos::RCP<Epetra_FECrsMatrix> GLpYUEpetraDataPool::getB() { return B_; }
500
501Teuchos::RCP<Epetra_FECrsMatrix> GLpYUEpetraDataPool::getH() { return H_; }
502
503Teuchos::RCP<Epetra_FECrsMatrix> GLpYUEpetraDataPool::getR() { return R_; }
504
505Teuchos::RCP<Epetra_CrsMatrix> GLpYUEpetraDataPool::getAugmat() { return Augmat_; }
506
507Teuchos::RCP<Epetra_FECrsMatrix> GLpYUEpetraDataPool::getNpy() { return Npy_; }
508
509Teuchos::RCP<Epetra_FEVector> GLpYUEpetraDataPool::getb() { return b_; }
510
511Teuchos::RCP<Epetra_FEVector> GLpYUEpetraDataPool::getq() { return q_; }
512
513Teuchos::RCP<Epetra_FEVector> GLpYUEpetraDataPool::getNy() { return Ny_; }
514
516
517Teuchos::RCP<const Epetra_SerialDenseMatrix> GLpYUEpetraDataPool::getipcoords() { return ipcoords_; }
518
519Teuchos::RCP<const Epetra_IntSerialDenseVector> GLpYUEpetraDataPool::getipindx() { return ipindx_; }
520
521Teuchos::RCP<const Epetra_SerialDenseMatrix> GLpYUEpetraDataPool::getpcoords() { return pcoords_; }
522
523Teuchos::RCP<const Epetra_IntSerialDenseVector> GLpYUEpetraDataPool::getpindx() { return pindx_; }
524
525Teuchos::RCP<const Epetra_IntSerialDenseMatrix> GLpYUEpetraDataPool::gett() { return t_; }
526
527Teuchos::RCP<const Epetra_IntSerialDenseMatrix> GLpYUEpetraDataPool::gete() { return e_; }
528
529
530void GLpYUEpetraDataPool::computeNy( const Teuchos::RCP<const Epetra_MultiVector> & y )
531{
532 Epetra_Map overlapmap(-1, pindx_->M(), (int*)(pindx_)->A(), 1, *commptr_);
533 Epetra_Map standardmap(A_->DomainMap());
534 Teuchos::RCP<Epetra_MultiVector> yoverlap = Teuchos::rcp(new Epetra_MultiVector(overlapmap, 1));
535 Epetra_Import Importer(overlapmap, standardmap);
536 yoverlap->Import(*y, Importer, Insert);
537 nonlinvec(*commptr_, *ipindx_, *ipcoords_, *pindx_, *pcoords_, *t_, yoverlap, Ny_);
538}
539
540
541void GLpYUEpetraDataPool::computeNpy( const Teuchos::RCP<const Epetra_MultiVector> & y )
542{
543 Epetra_Map overlapmap(-1, pindx_->M(), (int*)(pindx_)->A(), 1, *commptr_);
544 Epetra_Map standardmap(A_->DomainMap());
545 Teuchos::RCP<Epetra_MultiVector> yoverlap = Teuchos::rcp(new Epetra_MultiVector(overlapmap, 1));
546 Epetra_Import Importer(overlapmap, standardmap);
547 yoverlap->Import(*y, Importer, Insert);
548 nonlinjac(*commptr_, *ipindx_, *ipcoords_, *pindx_, *pcoords_, *t_, yoverlap, Npy_);
549}
550
551
553{
554 Epetra_Map standardmap(A_->DomainMap());
555 Epetra_Map bdryctrlmap(B_->DomainMap());
556
557 int indexBase = 1;
558
559 int numstates = standardmap.NumGlobalElements();
560 //int numcontrols = bdryctrlmap.NumGlobalElements();
561 int nummystates = standardmap.NumMyElements();
562 int nummycontrols = bdryctrlmap.NumMyElements();
563
564 Epetra_IntSerialDenseVector KKTmapindx(2*nummystates+nummycontrols);
565
566 // Build KKT map.
567 Epetra_IntSerialDenseVector states(nummystates);
568 Epetra_IntSerialDenseVector controls(nummycontrols);
569 standardmap.MyGlobalElements(states.Values());
570 bdryctrlmap.MyGlobalElements(controls.Values());
571 for (int i=0; i<nummystates; i++) {
572 KKTmapindx(i) = states(i);
573 KKTmapindx(nummystates+nummycontrols+i) = 2*numstates+states(i);
574 }
575 for (int i=0; i<nummycontrols; i++) {
576 KKTmapindx(nummystates+i) = numstates+controls(i);
577 }
578 Epetra_Map KKTmap(-1, 2*nummystates+nummycontrols, KKTmapindx.Values(), indexBase, *(commptr_));
579
580
581 // Start building the KKT matrix.
582
583 Augmat_ = Teuchos::rcp(new Epetra_CrsMatrix(Copy, KKTmap, 0));
584
585 double one[1];
586 one[0]=1.0;
587 for (int i=0; i<nummystates+nummycontrols; i++) {
588 Augmat_->InsertGlobalValues(KKTmapindx.Values()[i], 1, one, KKTmapindx.Values()+i);
589 }
590
591 int checkentries=0;
592 int nummyentries=0;
593 Epetra_SerialDenseVector values(nummyentries);
594 Epetra_IntSerialDenseVector indices(nummyentries);
595 // Insert A and Npy into Augmat.
596 for (int i=0; i<nummystates; i++) {
597 nummyentries = A_->NumMyEntries(i);
598 values.Resize(nummyentries);
599 indices.Resize(nummyentries);
600 A_->ExtractGlobalRowCopy(KKTmapindx.Values()[i], nummyentries, checkentries, values.Values(),
601 indices.Values());
602 if (nummyentries > 0)
603 Augmat_->InsertGlobalValues(KKTmapindx.Values()[i]+2*numstates, nummyentries, values.Values(),
604 indices.Values());
605 nummyentries = Npy_->NumMyEntries(i);
606 values.Resize(nummyentries);
607 indices.Resize(nummyentries);
608 Npy_->ExtractGlobalRowCopy(KKTmapindx.Values()[i], nummyentries, checkentries, values.Values(),
609 indices.Values());
610 if (nummyentries > 0)
611 Augmat_->InsertGlobalValues(KKTmapindx.Values()[i]+2*numstates, nummyentries, values.Values(),
612 indices.Values());
613 }
614 // Insert B into Augmat.
615 for (int i=0; i<nummystates; i++) {
616 nummyentries = B_->NumMyEntries(i);
617 values.Resize(nummyentries);
618 indices.Resize(nummyentries);
619 B_->ExtractGlobalRowCopy(KKTmapindx.Values()[i], nummyentries, checkentries, values.Values(),
620 indices.Values());
621 for (int j=0; j<nummyentries; j++)
622 indices[j] = indices[j]+numstates;
623 if (nummyentries > 0)
624 Augmat_->InsertGlobalValues(KKTmapindx.Values()[i]+2*numstates, nummyentries, values.Values(),
625 indices.Values());
626 }
627
629 Epetra_CrsMatrix & transA = dynamic_cast<Epetra_CrsMatrix&>(transposer(*A_));
630 Epetra_CrsMatrix & transB = dynamic_cast<Epetra_CrsMatrix&>(transposer(*B_));
631 Epetra_CrsMatrix & transNpy = dynamic_cast<Epetra_CrsMatrix&>(transposer(*Npy_));
632 // Insert transpose of A and Npy into Augmat.
633 for (int i=0; i<nummystates; i++) {
634 nummyentries = transA.NumMyEntries(i);
635 values.Resize(nummyentries);
636 indices.Resize(nummyentries);
637 transA.ExtractGlobalRowCopy(KKTmapindx.Values()[i], nummyentries, checkentries, values.Values(),
638 indices.Values());
639 for (int j=0; j<nummyentries; j++)
640 indices[j] = indices[j]+2*numstates;
641 if (nummyentries > 0)
642 Augmat_->InsertGlobalValues(KKTmapindx.Values()[i], nummyentries, values.Values(),
643 indices.Values());
644 nummyentries = transNpy.NumMyEntries(i);
645 values.Resize(nummyentries);
646 indices.Resize(nummyentries);
647 transNpy.ExtractGlobalRowCopy(KKTmapindx.Values()[i], nummyentries, checkentries, values.Values(),
648 indices.Values());
649 for (int j=0; j<nummyentries; j++)
650 indices[j] = indices[j]+2*numstates;
651 if (nummyentries > 0)
652 Augmat_->InsertGlobalValues(KKTmapindx.Values()[i], nummyentries, values.Values(),
653 indices.Values());
654 }
655 // Insert transpose of B into Augmat.
656 for (int i=0; i<nummystates; i++) {
657 nummyentries = transB.NumMyEntries(i);
658 values.Resize(nummyentries);
659 indices.Resize(nummyentries);
660 transB.ExtractGlobalRowCopy(KKTmapindx.Values()[i], nummyentries, checkentries, values.Values(),
661 indices.Values());
662 for (int j=0; j<nummyentries; j++)
663 indices[j] = indices[j]+2*numstates;
664 int err = 0;
665 if (nummyentries > 0)
666 err = Augmat_->InsertGlobalValues(KKTmapindx.Values()[i]+numstates, nummyentries,
667 values.Values(), indices.Values());
668 // This will give a nasty message if something goes wrong with the insertion of B transpose.
669 if (err < 0) {
670 std::cout << "Insertion of entries failed:\n";
671 std::cout << indices;
672 std::cout << nummyentries << std::endl;
673 std::cout << "at row: " << KKTmapindx.Values()[i]+numstates << std::endl << std::endl;
674 }
675 }
676
677 Augmat_->FillComplete(KKTmap, KKTmap);
678 // End building the KKT matrix.
679
680}
681
682void GLpYUEpetraDataPool::PrintVec( const Teuchos::RCP<const Epetra_Vector> & x )
683{
684 Vector2MATLAB(*x, std::cout);
685}
686
688 Epetra_BLAS blas;
689 Epetra_SerialDenseVector product(4);
690
691 // get nodes and weights
693
694 // Evaluate nodal basis functions and their derivatives at quadrature
695 // pts N(i,j) = value of the j-th basis function at quadrature node i.
696 N.Shape(Nodes.M(),3);
697 for (int i=0; i<Nodes.M(); i++) {
698 N(i,0) = 1.0 - Nodes(i,0) - Nodes(i,1);
699 N(i,1) = Nodes(i,0);
700 N(i,2) = Nodes(i,1);
701 }
702 // Nx1(i,j) partial derrivative of basis function j wrt x1 at node i
703 Nx1.Shape(Nodes.M(),3);
704 for (int i=0; i<Nodes.M(); i++) {
705 Nx1(i,0) = -1.0;
706 Nx1(i,1) = 1.0;
707 Nx1(i,2) = 0.0;
708 }
709 // Nx2(i,j) partial derrivative of basis function j wrt x2 at node i
710 Nx2.Shape(Nodes.M(),3);
711 for (int i=0; i<Nodes.M(); i++) {
712 Nx2(i,0) = -1.0;
713 Nx2(i,1) = 0.0;
714 Nx2(i,2) = 1.0;
715 }
716
717 // Evaluate integrals of various combinations of partial derivatives
718 // of the local basis functions (they're constant).
719 S1.Shape(3,3);
720 S1(0,0)= 1.0; S1(0,1)=-1.0; S1(0,2)= 0.0;
721 S1(1,0)=-1.0; S1(1,1)= 1.0; S1(1,2)= 0.0;
722 S1(2,0)= 0.0; S1(2,1)= 0.0; S1(2,2)= 0.0;
723 S2.Shape(3,3);
724 S2(0,0)= 2.0; S2(0,1)=-1.0; S2(0,2)=-1.0;
725 S2(1,0)=-1.0; S2(1,1)= 0.0; S2(1,2)= 1.0;
726 S2(2,0)=-1.0; S2(2,1)= 1.0; S2(2,2)= 0.0;
727 S3.Shape(3,3);
728 S3(0,0)= 1.0; S3(0,1)= 0.0; S3(0,2)=-1.0;
729 S3(1,0)= 0.0; S3(1,1)= 0.0; S3(1,2)= 0.0;
730 S3(2,0)=-1.0; S3(2,1)= 0.0; S3(2,2)= 1.0;
731
732 // Evaluate integrals of basis functions N(i).
733 Nw.Size(3);
734 Nw.Multiply('T', 'N', 1.0, N, Weights, 0.0);
735
736 // Evaluate integrals of 9 products N(i)*N(j).
737 NNw.Shape(3,3);
738 for (int i=0; i<3; i++) {
739 for (int j=0; j<3; j++) {
740 compproduct(product, N[i], N[j]);
741 NNw(i,j) = blas.DOT(Weights.M(), Weights.A(), product.A());
742 }
743 }
744
745 // Evaluate integrals of 27 products N(i)*N(j)*N(k).
747 NNNw[0].Shape(3,3); NNNw[1].Shape(3,3); NNNw[2].Shape(3,3);
748 for (int i=0; i<3; i++) {
749 for (int j=0; j<3; j++) {
750 for (int k=0; k<3; k++) {
751 compproduct(product, N[i], N[j], N[k]);
752 NNNw[k](i,j) = blas.DOT(Weights.M(), Weights.A(), product.A());
753 }
754 }
755 }
756
757 // Evaluate integrals of 27 products N(i)*(dN(j)/dx1)*N(k).
759 NdNdx1Nw[0].Shape(3,3); NdNdx1Nw[1].Shape(3,3); NdNdx1Nw[2].Shape(3,3);
760 for (int i=0; i<3; i++) {
761 for (int j=0; j<3; j++) {
762 for (int k=0; k<3; k++) {
763 compproduct(product, N[i], Nx1[j], N[k]);
764 NdNdx1Nw[k](i,j) = blas.DOT(Weights.M(), Weights.A(), product.A());
765 }
766 }
767 }
768
769 // Evaluate integrals of 27 products N(i)*(dN(j)/dx2)*N(k).
771 NdNdx2Nw[0].Shape(3,3); NdNdx2Nw[1].Shape(3,3); NdNdx2Nw[2].Shape(3,3);
772 for (int i=0; i<3; i++) {
773 for (int j=0; j<3; j++) {
774 for (int k=0; k<3; k++) {
775 compproduct(product, N[i], Nx2[j], N[k]);
776 NdNdx2Nw[k](i,j) = blas.DOT(Weights.M(), Weights.A(), product.A());
777 }
778 }
779 }
780
781}
782
783void Usr_Par::Print(std::ostream& os) const {
784 os << std::endl;
785 os << "\n\nQuadrature nodes:";
786 os << "\n-----------------";
787 Nodes.Print(os);
788 os << "\n\nQuadrature weights:";
789 os << "\n-------------------\n";
790 Weights.Print(os);
791 os << "\n\nNodal basis functions:";
792 os << "\n----------------------";
793 N.Print(os);
794 os << "\n\nIntegrals of combinations of partial derivatives:";
795 os << "\n-------------------------------------------------";
796 S1.Print(os);
797 S2.Print(os);
798 S3.Print(os);
799 os << "\n\nIntegrals of basis functions:";
800 os << "\n-----------------------------\n";
801 Nw.Print(os);
802 os << "\n\nIntegrals of products N(i)*N(j):";
803 os << "\n--------------------------------\n";
804 NNw.Print(os);
805 os << "\n\nIntegrals of products N(i)*N(j)*N(k):";
806 os << "\n-------------------------------------\n";
807 NNNw[0].Print(os); NNNw[1].Print(os); NNNw[2].Print(os);
808 os << "\n\nIntegrals of products N(i)*(dN(j)/dx1)*N(k):";
809 os << "\n--------------------------------------------\n";
810 NdNdx1Nw[0].Print(os); NdNdx1Nw[1].Print(os); NdNdx1Nw[2].Print(os);
811 os << "\n\nIntegrals of products N(i)*(dN(j)/dx2)*N(k):";
812 os << "\n--------------------------------------------\n";
813 NdNdx2Nw[0].Print(os); NdNdx2Nw[1].Print(os); NdNdx2Nw[2].Print(os);
814}
815
816std::ostream& operator <<(std::ostream & out, const Usr_Par & usr_par)
817{
818 usr_par.Print(out);
819 return out;
820}
821
822} // namespace GLpApp
823
824//
825// Implementations
826//
827
829 double *first, double *second)
830{
831 for (int i=0; i<product.M(); i++) {
832 product[i] = first[i]*second[i];
833 }
834 return(0);
835}
836
838 double *first, double *second, double *third)
839{
840 for (int i=0; i<product.M(); i++) {
841 product[i] = first[i]*second[i]*third[i];
842 }
843 return(0);
844}
845
846//#define GLPAPP_SHOW_BOUNDARY_ASSEMBLY
847
848/* \brief Performs finite-element assembly of face mass matrices.
849
850 \param Comm [in] - The Epetra (MPI) communicator.
851 \param ipindx [in] - Vector of NUMIP indices of nodes that are `unique' to a subdomain
852 (i.e. owned by the corresponding processor).
853 \param ipcoords [in] - Matrix (NUMIP x 2) of x-y-coordinates of the vertices cooresponding
854 to indices ipindx: \n
855 ipcoords(i,0) x-coordinate of node i, \n
856 ipcoords(i,1) y-coordinate of node i.
857 \param pindx [in] - Vector of NUMP indices of ALL nodes in a subdomain, including
858 the shared nodes.
859 \param pcoords [in] - Matrix (NUMP x 2) of x-y-coordinates of the vertices cooresponding
860 to indices pindx: \n
861 pcoords(i,0) x-coordinate of node i, \n
862 pcoords(i,1) y-coordinate of node i.
863 \param t [in] - Matrix (ELE x 3) of indices of the vertices in a triangle: \n
864 t(i,j) index of the 0-based j-th vertex in triangle i, where i = 0, ..., numElements-1
865 \param e [in] - Matrix (EDGE x 3) of edges. \n
866 e(i,1:2) contains the indices of the endpoints of edge i, where i = 0, ..., numEdges-1 \n
867 e(i,3) contains the boundary marker
868 \param B_out [out] - Reference-counting pointer to the Epetra_FECrsMatrix describing the FE
869 state/control face mass matrix.
870 \param R_out [out] - Reference-counting pointer to the Epetra_FECrsMatrix describing the FE
871 control/control volume mass matrix.
872 \return 0 if successful.
873
874 \par Detailed Description:
875
876 -# Assembles the FE state/control mass matrix \e B, given by
877 \f[
878 \mathbf{B}_{jk} = b(\mu_k,\phi_j) = - \int_{\partial \Omega} \mu_k(x) \phi_j(x) dx,
879 \f]
880 where \f$\{ \phi_j \}_{j = 1}^{m}\f$ is the piecewise linear nodal basis for the finite-dimensional
881 state space, and \f$\{ \mu_j \}_{j = 1}^{n}\f$ is the piecewise linear nodal basis for the
882 finite-dimensional control space.
883 -# Assembles the FE control/control mass matrix \e R, given by
884 \f[
885 \mathbf{R}_{jk} = b(\mu_k,\mu_j) = - \int_{\partial \Omega} \mu_k(x) \mu_j(x) dx,
886 \f]
887 where \f$\{ \mu_j \}_{j = 1}^{n}\f$ is the piecewise linear nodal basis for the finite-dimensional
888 control space.
889*/
891 const Epetra_Comm &Comm
892 ,const Epetra_IntSerialDenseVector &ipindx
893 ,const Epetra_SerialDenseMatrix &ipcoords
894 ,const Epetra_IntSerialDenseVector &pindx
895 ,const Epetra_SerialDenseMatrix &pcoords
898 ,Teuchos::RCP<Epetra_FECrsMatrix> *B_out
899 ,Teuchos::RCP<Epetra_FECrsMatrix> *R_out
900 )
901{
902
903 using Teuchos::rcp;
904
905#ifdef GLPAPP_SHOW_BOUNDARY_ASSEMBLY
906 Teuchos::RCP<Teuchos::FancyOStream>
907 out = Teuchos::VerboseObjectBase::getDefaultOStream();
908 Teuchos::OSTab tab(out);
909 *out << "\nEntering assemble_bdry(...) ...\n";
910#endif
911
912 int numLocElems = t.M();
913 int numLocEdges = e.M();
914
915 int indexBase = 1;
916
917 //
918 // Create a sorted (by global ID) list of boundry nodes
919 //
920 int * lastindx = 0;
921 Epetra_IntSerialDenseVector BdryNode(2*numLocEdges);
922 for (int j=0; j<numLocEdges; j++) {
923 BdryNode[j] = e(j,0);
924 BdryNode[j+numLocEdges] = e(j,1);
925 }
926 std::sort(BdryNode.Values(), BdryNode.Values()+2*numLocEdges);
927 lastindx = std::unique(BdryNode.Values(), BdryNode.Values()+2*numLocEdges);
928 const int numBdryNodes = lastindx - BdryNode.Values();
929 BdryNode.Resize(numBdryNodes); // RAB: This does not overwrite?
930
931 //
932 // Determine the boundary vertices that belong to this processor.
933 //
934 Epetra_IntSerialDenseVector MyBdryNode(numBdryNodes);
935 lastindx = std::set_intersection(
936 BdryNode.Values(), BdryNode.Values()+numBdryNodes, // (Sorted) set A
937 ipindx.Values(), ipindx.Values()+ipindx.M(), // (Sorted) set B
938 MyBdryNode.Values() // Intersection
939 );
940 const int numMyBdryNodes = lastindx - MyBdryNode.Values();
941 MyBdryNode.Resize(numMyBdryNodes); // RAB: This does not overwrite?
942
943 //
944 // Define the maps for the various lists
945 //
946 Epetra_Map standardmap(-1, ipindx.M(), const_cast<int*>(ipindx.A()), indexBase, Comm);
947 Epetra_Map overlapmap(-1, pindx.M(), const_cast<int*>(pindx.A()), indexBase, Comm);
948 Epetra_Map mybdryctrlmap(-1, numMyBdryNodes, const_cast<int*>(MyBdryNode.A()), indexBase, Comm);
949 // Above, it is important to note what mybndyctrlmap represents. It is the
950 // a sorted list of global vertex node IDS for nodes on a boundary that are
951 // uniquely owned by the local process.
952
953#ifdef GLPAPP_SHOW_BOUNDARY_ASSEMBLY
954 *out << "\nstandardmap:\n";
955 standardmap.Print(Teuchos::OSTab(out).o());
956 *out << "\nmybdryctrlmap:\n";
957 mybdryctrlmap.Print(Teuchos::OSTab(out).o());
958#endif
959
960 //
961 // Allocate matrices to fill
962 //
963 Teuchos::RCP<Epetra_FECrsMatrix>
964 B = rcp(new Epetra_FECrsMatrix(Copy,standardmap,0)),
965 R = rcp(new Epetra_FECrsMatrix(Copy,standardmap,0));
966 // NOTE: The data map is the same as for the volume matrices. Later, when
967 // FillComplete is called, we will fix the proper domain and range maps.
968 // Declare quantities needed for the call to the local assembly routine.
969 const int numNodesPerEdge = 2;
970 Epetra_IntSerialDenseVector epetra_nodes(numNodesPerEdge);
971
972 //
973 // Load B and R by looping through the edges
974 //
975
977 int err=0;
978
979 for( int i=0; i < numLocEdges; i++ ) {
980
981 const int
982 global_id_0 = e(i,0),
983 global_id_1 = e(i,1),
984 local_id_0 = overlapmap.LID(global_id_0), // O(log(numip)) bindary search
985 local_id_1 = overlapmap.LID(global_id_1); // O(log(numip)) bindary search
986
987 epetra_nodes(0) = global_id_0;
988 epetra_nodes(1) = global_id_1;
989
990 const double
991 x0 = pcoords(local_id_0,0),
992 y0 = pcoords(local_id_0,1),
993 x1 = pcoords(local_id_1,0),
994 y1 = pcoords(local_id_1,1);
995
996 const double l = sqrt(pow(x0-x1,2) + pow(y0-y1,2)); // Length of this edge
997
998 // We have an explicit formula for Bt.
999 const double l_sixth = l * (1.0/6.0);
1000 Bt(0,0) = l_sixth * 2.0;
1001 Bt(0,1) = l_sixth * 1.0;
1002 Bt(1,0) = l_sixth * 1.0;
1003 Bt(1,1) = l_sixth * 2.0;
1004
1005#ifdef GLPAPP_SHOW_BOUNDARY_ASSEMBLY
1006 *out
1007 << "\nEdge global nodes = ("<<global_id_0<<","<<global_id_1<<"),"
1008 << " local nodes = ("<<local_id_0<<","<<local_id_1<<"),"
1009 << " (x0,y0)(x1,y1)=("<<x0<<","<<y0<<")("<<x1<<","<<y1<<"),"
1010 << " Bt = ["<<Bt(0,0)<<","<<Bt(0,1)<<";"<<Bt(1,0)<<","<<Bt(1,1)<<"]\n";
1011#endif
1012
1013 const int format = Epetra_FECrsMatrix::COLUMN_MAJOR;
1014 err = B->InsertGlobalValues(epetra_nodes,Bt,format);
1015 if (err<0) return(err);
1016 err = R->InsertGlobalValues(epetra_nodes,Bt,format);
1017 if (err<0) return(err);
1018
1019 }
1020
1021/*
1022
1023 err = B->GlobalAssemble(false);
1024 if (err<0) return(err);
1025 err = R->GlobalAssemble(false);
1026 if (err<0) return(err);
1027
1028 err = B->FillComplete(mybdryctrlmap,standardmap);
1029 if (err<0) return(err);
1030 err = R->FillComplete(mybdryctrlmap,mybdryctrlmap);
1031 if (err<0) return(err);
1032
1033*/
1034
1035 err = B->GlobalAssemble(mybdryctrlmap,standardmap,true);
1036 if (err<0) return(err);
1037 err = R->GlobalAssemble(mybdryctrlmap,mybdryctrlmap,true);
1038 if (err<0) return(err);
1039
1040 if(B_out) *B_out = B;
1041 if(R_out) *R_out = R;
1042
1043#ifdef GLPAPP_SHOW_BOUNDARY_ASSEMBLY
1044 *out << "\nB =\n";
1045 B->Print(Teuchos::OSTab(out).o());
1046 *out << "\nLeaving assemble_bdry(...) ...\n";
1047#endif
1048
1049 return(0);
1050
1051}
1052
1053/* \brief Performs finite-element assembly of volume stiffness and mass matrices,
1054 and the right-hand side (forcing term).
1055
1056 \param Comm [in] - The Epetra (MPI) communicator.
1057 \param ipindx [in] - Vector of NUMIP indices of nodes that are `unique' to a subdomain
1058 (i.e. owned by the corresponding processor).
1059 \param ipcoords [in] - Matrix (NUMIP x 2) of x-y-coordinates of the vertices cooresponding
1060 to indices ipindx: \n
1061 ipcoords(i,0) x-coordinate of node i, \n
1062 ipcoords(i,1) y-coordinate of node i.
1063 \param pindx [in] - Vector of NUMP indices of ALL nodes in a subdomain, including
1064 the shared nodes.
1065 \param pcoords [in] - Matrix (NUMP x 2) of x-y-coordinates of the vertices cooresponding
1066 to indices pindx: \n
1067 pcoords(i,0) x-coordinate of node i, \n
1068 pcoords(i,1) y-coordinate of node i.
1069 \param t [in] - Matrix (ELE x 3) of indices of the vertices in a triangle: \n
1070 t(i,j) index of the j-th vertex in triangle i, where i = 1, ..., ELE
1071 \param e [in] - Matrix (EDGE x 3) of edges. \n
1072 e(i,1:2) contains the indices of the endpoints of edge i, where i = 1, ..., EDGE \n
1073 e(i,3) contains the boundary marker \n
1074 e(i,3) = 1 Dirichlet boundary conditions on edge i \n
1075 e(i,3) = 2 Neumann boundary conditions on edge i \n
1076 e(i,3) = 3 Robin boundary conditions on edge i
1077 \param A [out] - Reference-counting pointer to the Epetra_FECrsMatrix describing the FE volume
1078 stiffness matrix for the advection-diffusion equation. Includes advection,
1079 diffusion, and reaction terms, and modifications that account for the boundary
1080 conditions.
1081 \param M [out] - Reference-counting pointer to the Epetra_FECrsMatrix describing the FE volume
1082 mass matrix.
1083 \param b [out] - Reference-counting pointer to the Epetra_FEVector describing the discretized
1084 forcing term in the advection-diffusion equation. Includes modifications that
1085 account for the boundary conditions.
1086 \return 0 if successful.
1087
1088 \par Detailed Description:
1089
1090 -# Assembles the FE volume stiffness matrix \e A and the right-hand side \e b for the
1091 solution of an advection-diffusion equation using piecewise linear finite elements.
1092 The advection-diffusion equation is
1093 \f{align*}
1094 - \nabla \cdot \left( K(x) \nabla y(x) \right) + \mathbf{c}(x) \cdot \nabla y(x) + r(x) y(x)
1095 &= f(x), & x &\in \Omega,\; \\
1096 y(x) &= d(x), & x &\in {\partial \Omega}_d, \\
1097 K(x) \frac{\partial}{\partial \mathbf{n}} y(x) &= g(x), & x &\in {\partial \Omega}_n, \\
1098 \sigma_0 K(x) \frac{\partial}{\partial \mathbf{n}} y(x)
1099 + \sigma_1 y(x) &= g(x), & x &\in {\partial \Omega}_r,
1100 \f}
1101 where \f$ K \f$ represents scalar diffusion, \f$ \mathbf{c} \f$ is the advection vector field,
1102 \f$ r \f$ is the reaction term, \f$ d \f$ and \f$ g \f$ are given functions, \f$sigma_0\f$ and
1103 \f$ sigma_1 \f$ are given quantities, \f$ {\partial \Omega}_d \f$ is the Dirichlet boundary,
1104 \f$ {\partial \Omega}_n \f$ is the Neumann boundary, and \f$ {\partial \Omega}_r \f$ is the Robin
1105 boundary. The quantities \f$ K \f$, \f$ \mathbf{c} \f$, \f$ r \f$, \f$ d \f$, and \f$ g \f$ are
1106 assumed to be piecewise linear. Currently, they are to be hard-coded inside this function.
1107 -# Assembles the FE volume mass matrix \e M.
1108*/
1109int GLpApp::assemble(const Epetra_Comm & Comm,
1110 const Epetra_IntSerialDenseVector & ipindx,
1111 const Epetra_SerialDenseMatrix & ipcoords,
1112 const Epetra_IntSerialDenseVector & pindx,
1113 const Epetra_SerialDenseMatrix & pcoords,
1116 RCP<Epetra_FECrsMatrix> & A,
1117 RCP<Epetra_FECrsMatrix> & M,
1118 RCP<Epetra_FEVector> & b)
1119{
1120
1121 int myPID = Comm.MyPID();
1122 int numProcs = Comm.NumProc();
1123 Usr_Par usr_par;
1124
1125 int numLocElems = t.M();
1126 int numNodesPerElem = 3;
1127
1128 int indexBase = 1;
1129
1130 Epetra_Map standardmap(-1, ipindx.M(), (int*)ipindx.A(), indexBase, Comm);
1131 Epetra_Map overlapmap(-1, pindx.M(), (int*)pindx.A(), indexBase, Comm);
1132
1133 int* nodes = new int[numNodesPerElem];
1134 int i=0, j=0, err=0;
1135
1136 A = rcp(new Epetra_FECrsMatrix(Copy, standardmap, 0));
1137 M = rcp(new Epetra_FECrsMatrix(Copy, standardmap, 0));
1138 b = rcp(new Epetra_FEVector(standardmap,1));
1139
1140 // Declare quantities needed for the call to the local assembly routine.
1142 Epetra_IntSerialDenseVector epetra_nodes(View, nodes, numNodesPerElem);
1143
1144
1145 /* ************************ First, we build A and b. ************************ */
1148 Epetra_SerialDenseMatrix vertices(numNodesPerElem, pcoords.N());
1149
1150 Epetra_SerialDenseVector k(numNodesPerElem);
1151 for (i=0; i< numNodesPerElem; i++) k(i)=1.0;
1152 Epetra_SerialDenseMatrix c(numNodesPerElem,2);
1153 for (i=0; i< numNodesPerElem; i++) { c(i,0)=0.0; c(i,1)=0.0; }
1154 Epetra_SerialDenseVector r(numNodesPerElem);
1155 for (i=0; i< numNodesPerElem; i++) r(i)=0.0;
1156 Epetra_SerialDenseVector f(numNodesPerElem);
1157 for (i=0; i< numNodesPerElem; i++) f(i)=0.0;
1159 g(0)=0.0; g(1)=0.0;
1160 Epetra_SerialDenseVector sigma(2);
1161 sigma(0)=0.0; sigma(1)=0.0;
1162 for(i=0; i<numLocElems; i++) {
1163 nodes[0] = t(i,0); nodes[1] = t(i,1); nodes[2] = t(i,2);
1164 for (j=0; j<numNodesPerElem; j++) {
1165 // get vertex
1166 vertices(j,0) = pcoords(overlapmap.LID(nodes[j]), 0);
1167 vertices(j,1) = pcoords(overlapmap.LID(nodes[j]), 1);
1168 // set rhs function
1169 f(j) = cos(GLp_pi*vertices(j,0))*cos(GLp_pi*vertices(j,1)) *
1170 (2*GLp_pi*GLp_pi + pow(cos(GLp_pi*vertices(j,0)),2)*pow(cos(GLp_pi*vertices(j,1)),2) - 1.0);
1171 }
1172 lassembly(vertices, k, c, r, f, usr_par, At, bt);
1173 err = A->InsertGlobalValues(epetra_nodes, At, format);
1174 if (err<0) return(err);
1175 err = b->SumIntoGlobalValues(epetra_nodes, bt);
1176 if (err<0) return(err);
1177 }
1178
1179 // MAKE ADJUSTMENTS TO A and b TO ACCOUNT FOR BOUNDARY CONDITIONS.
1180
1181 // Get Neumann boundary edges.
1182 Epetra_IntSerialDenseMatrix neumann(e.M(),2);
1183 j = 0;
1184 for (i=0; i<e.M(); i++) {
1185 if (e(i,2)==2) {
1186 neumann(j,0) = e(i,0); neumann(j,1) = e(i,1);
1187 j++;
1188 }
1189 }
1190 neumann.Reshape(j,2);
1191 // Adjust for Neumann BC's.
1192 if (neumann.M() != 0) {
1193 Epetra_BLAS blas;
1194 Epetra_SerialDenseMatrix quadnodes;
1195 Epetra_SerialDenseVector quadweights;
1198 Epetra_SerialDenseVector product(2);
1199
1200 // Get quadrature weights and points.
1201 quadrature(1, 2, quadnodes, quadweights);
1202
1203 // Evaluate nodal basis functions at quadrature points
1204 // N(i,j) value of basis function j at quadrature node i
1205 N.Shape(quadnodes.M(),2);
1206 for (i=0; i<quadnodes.M(); i++) {
1207 N(i,0) = 1.0 - quadnodes(i,0);
1208 N(i,1) = quadnodes(i,0);
1209 }
1210
1211 // Evaluate integrals of 4 products N(i)*N(j).
1212 NN.Shape(2,2);
1213 for (i=0; i<2; i++) {
1214 for (j=0; j<2; j++) {
1215 compproduct(product, N[i], N[j]);
1216 NN(i,j) = blas.DOT(quadweights.M(), quadweights.A(), product.A());
1217 }
1218 }
1219
1220 Epetra_IntSerialDenseVector neumann_nodes(2);
1221 Epetra_SerialDenseVector neumann_add(2);
1222 double l;
1223 for (i=0; i<neumann.M(); i++) {
1224 neumann_nodes(0) = neumann(i,0); neumann_nodes(1) = neumann(i,1);
1225 neumann_add(0) = pcoords(overlapmap.LID(neumann_nodes(0)),0)
1226 -pcoords(overlapmap.LID(neumann_nodes(1)),0);
1227 neumann_add(1) = pcoords(overlapmap.LID(neumann_nodes(0)),1)
1228 -pcoords(overlapmap.LID(neumann_nodes(1)),1);
1229 l = blas.NRM2(neumann_add.M(), neumann_add.A());
1230 neumann_add.Multiply('N', 'N', 1.0, NN, g, 0.0);
1231 neumann_add.Scale(l);
1232 err = b->SumIntoGlobalValues(neumann_nodes, neumann_add);
1233 if (err<0) return(err);
1234 }
1235 }
1236
1237 // Get Robin boundary edges.
1238 Epetra_IntSerialDenseMatrix robin(e.M(),2);
1239 j = 0;
1240 for (i=0; i<e.M(); i++) {
1241 if (e(i,2)==3) {
1242 robin(j,0) = e(i,0); robin(j,1) = e(i,1);
1243 j++;
1244 }
1245 }
1246 robin.Reshape(j,2);
1247 // Adjust for Robin BC's.
1248 if (robin.M() != 0) {
1249 Epetra_BLAS blas;
1250 Epetra_SerialDenseMatrix quadnodes;
1251 Epetra_SerialDenseVector quadweights;
1255 Epetra_SerialDenseVector product(2);
1256
1257 // Get quadrature weights and points.
1258 quadrature(1, 2, quadnodes, quadweights);
1259
1260 // Evaluate nodal basis functions at quadrature points
1261 // N(i,j) value of basis function j at quadrature node i
1262 N.Shape(quadnodes.M(),2);
1263 for (i=0; i<quadnodes.M(); i++) {
1264 N(i,0) = 1.0 - quadnodes(i,0);
1265 N(i,1) = quadnodes(i,0);
1266 }
1267
1268 // Evaluate integrals of 4 products N(i)*N(j).
1269 NN.Shape(2,2);
1270 for (i=0; i<2; i++) {
1271 for (j=0; j<2; j++) {
1272 compproduct(product, N[i], N[j]);
1273 NN(i,j) = blas.DOT(quadweights.M(), quadweights.A(), product.A());
1274 }
1275 }
1276
1277 // Evaluate integrals of 8 products N(i)*N(j)*N(k).
1278 NNN = new Epetra_SerialDenseMatrix[2];
1279 NNN[0].Shape(2,2); NNN[1].Shape(2,2);
1280 for (i=0; i<2; i++) {
1281 for (j=0; j<2; j++) {
1282 for (int k=0; k<2; k++) {
1283 compproduct(product, N[i], N[j], N[k]);
1284 NNN[k](i,j) = blas.DOT(quadweights.M(), quadweights.A(),
1285 product.A());
1286 }
1287 }
1288 }
1289
1290 Epetra_IntSerialDenseVector robin_nodes(2);
1291 Epetra_SerialDenseVector robin_b_add(2);
1292 Epetra_SerialDenseMatrix robin_A_add(2,2);
1293 double l;
1294 for (i=0; i<robin.M(); i++) {
1295 robin_nodes(0) = robin(i,0); robin_nodes(1) = robin(i,1);
1296
1297 robin_b_add(0) = pcoords(overlapmap.LID(robin_nodes(0)),0)
1298 -pcoords(overlapmap.LID(robin_nodes(1)),0);
1299 robin_b_add(1) = pcoords(overlapmap.LID(robin_nodes(0)),1)
1300 -pcoords(overlapmap.LID(robin_nodes(1)),1);
1301 l = blas.NRM2(robin_b_add.M(), robin_b_add.A());
1302 robin_b_add.Multiply('N', 'N', 1.0, NN, g, 0.0);
1303 robin_b_add.Scale(l);
1304 err = b->SumIntoGlobalValues(robin_nodes, robin_b_add);
1305 if (err<0) return(err);
1306
1307 NNN[0].Scale(sigma(0)); NNN[1].Scale(sigma(1));
1308 robin_A_add += NNN[0]; robin_A_add += NNN[1];
1309 robin_A_add.Scale(l);
1310 err = A->InsertGlobalValues(robin_nodes, robin_A_add, format);
1311 if (err<0) return(err);
1312 }
1313
1314 delete [] NNN;
1315 }
1316
1317 // Get Dirichlet boundary edges.
1318 Epetra_IntSerialDenseMatrix dirichlet(e.M(),2);
1319 j = 0;
1320 for (i=0; i<e.M(); i++) {
1321 if (e(i,2)==1) {
1322 dirichlet(j,0) = e(i,0); dirichlet(j,1) = e(i,1);
1323 j++;
1324 }
1325 }
1326 dirichlet.Reshape(j,2);
1327 // DIRICHLET NOT DONE! DO THIS LATER!!!!
1328
1329 /* ************************ Done building A and b. ************************ */
1330
1331
1332
1333 /* ************************ Second, we build M. ************************ */
1334
1336
1337 for (i=0; i< numNodesPerElem; i++) k(i)=0.0;
1338 for (i=0; i< numNodesPerElem; i++) { c(i,0)=0.0; c(i,1)=0.0; }
1339 for (i=0; i< numNodesPerElem; i++) r(i)=1.0;
1340 for (i=0; i< numNodesPerElem; i++) f(i)=0.0;
1341 g(0)=0.0; g(1)=0.0;
1342 sigma(0)=0.0; sigma(1)=0.0;
1343 for(i=0; i<numLocElems; i++) {
1344 nodes[0] = t(i,0); nodes[1] = t(i,1); nodes[2] = t(i,2);
1345 for (j=0; j<numNodesPerElem; j++) {
1346 vertices(j,0) = pcoords(overlapmap.LID(nodes[j]), 0);
1347 vertices(j,1) = pcoords(overlapmap.LID(nodes[j]), 1);
1348 }
1349 lassembly(vertices, k, c, r, f, usr_par, Mt, bt);
1350 err = M->InsertGlobalValues(epetra_nodes, Mt, format);
1351 if (err<0) return(err);
1352 }
1353
1354 /* ************************ Done building M. ************************ */
1355
1356
1357
1358 // Call global assemble and FillComplete at the same time.
1359
1360 err = A->GlobalAssemble();
1361 if (err<0) return(err);
1362
1363 err = b->GlobalAssemble();
1364 if (err<0) return(err);
1365
1366 err = M->GlobalAssemble();
1367 if (err<0) return(err);
1368
1369 delete [] nodes;
1370
1371 return(0);
1372}
1373
1374/* \brief Computes local stiffness matrix and local RHS vector for simplex
1375 (triangles in two dimensions).
1376
1377 \param vertices [in] - Matrix containing the global coordinates of the current simplex.
1378 \param k [in] - Vector containing the value of the diffusion \f$k\f$ at each vertex
1379 (\f$k\f$ is assumed to be piecewise linear), where \f$k\f$ is the coeff in the
1380 term \f$ \nabla \cdot (k \nabla(u)) \f$ in the advection-diffusion equation.
1381 \param c [in] - Matrix containing the value of the advection field \f$ \mathbf{c} \f$ at each
1382 vertex (\f$ \mathbf{c} \f$ is assumed to be piecewise linear), where
1383 \f$ \mathbf{c} \f$ is the 2-vector of coeffs in the term
1384 \f$ \mathbf{c}(x) \cdot \nabla y(x) \f$ in the advection-diffusion equation.
1385 \param r [in] - Vector containing the value of \f$ r \f$ at each vertex (\f$ r \f$ is assumed
1386 to be piecewise linear), where \f$ r \f$ is the coefficient in the term
1387 \f$ ru \f$ in the advection-diffusion equation.
1388 \param f [in] - Vector containing the value of \f$ f \f$ at each vertex (\f$ f \f$ is assumed to be
1389 piecewise linear), where \f$ f \f$ is the right hand side in the adv-diff eq
1390 \param usr_par [in] - class containing: \n
1391 - S1, S2, S3 (3x3) the integrals of various combinations of partials
1392 of local basis functions
1393 - N (1x3) integrals of local basis functions
1394 - NNN[3] (3x3), integrals of products of local basis functions N(i)*N(j)*N(k)
1395 - etc.
1396 \param At [out] - stiffness matrix contribution for the simplex
1397 \param bt [out] - right-hand side contribution for the simplex
1398
1399 \return 0 if successful.
1400
1401 \par Detailed Description:
1402
1403 Computes the local (per simplex) contributions to the FE volume stiffness matrix \e A and the
1404 right-hand side \e b for the advection-diffusion equation
1405 \f{align*}
1406 - \nabla \cdot \left( K(x) \nabla y(x) \right) + \mathbf{c}(x) \cdot \nabla y(x) + r(x) y(x)
1407 &= f(x), & x &\in \Omega,\; \\
1408 y(x) &= d(x), & x &\in {\partial \Omega}_d, \\
1409 K(x) \frac{\partial}{\partial \mathbf{n}} y(x) &= g(x), & x &\in {\partial \Omega}_n, \\
1410 \sigma_0 K(x) \frac{\partial}{\partial \mathbf{n}} y(x)
1411 + \sigma_1 y(x) &= g(x), & x &\in {\partial \Omega}_r.
1412 \f}
1413*/
1415 const Epetra_SerialDenseVector & k,
1416 const Epetra_SerialDenseMatrix & c,
1417 const Epetra_SerialDenseVector & r,
1418 const Epetra_SerialDenseVector & f,
1419 const Usr_Par & usr_par,
1422{
1423 // Note that the constructors below initialize all entries to 0.
1426 Epetra_SerialDenseMatrix BtB(2,2);
1431 Epetra_SerialDenseMatrix tmp(3,3);
1432 double dB, adB;
1433 At.Shape(3,3);
1434 bt.Size(3);
1435
1436 // Construct affine transformation matrix.
1437 for(int i=0; i<2; i++) {
1438 B(i,0) = vertices(1,i)-vertices(0,i);
1439 B(i,1) = vertices(2,i)-vertices(0,i);
1440 }
1441 dB = determinant(B);
1442 adB = abs(dB);
1443
1444 // Compute matrix C = inv(B'*B).
1445 BtB.Multiply('T', 'N', 1.0, B, B, 0.0);
1446 inverse(BtB, C);
1447
1448 inverse(B, b); b.Scale(dB);
1449
1450 // Compute the part corresponding to div(K*grad(u)).
1451 tmp = usr_par.S1; tmp.Scale(C(0,0));
1452 M1 += tmp;
1453 tmp = usr_par.S2; tmp.Scale(C(0,1));
1454 M1 += tmp;
1455 tmp = usr_par.S3; tmp.Scale(C(1,1));
1456 M1 += tmp;
1457 M1.Scale( (k(0)*usr_par.Nw(0) + k(1)*usr_par.Nw(1) +
1458 k(2)*usr_par.Nw(2)) * adB );
1459
1460 // Compute the part corresponding to c'*grad(u).
1461 tmp = usr_par.NdNdx1Nw[0]; tmp.Scale(b(0,0)*c(0,0)+b(0,1)*c(0,1));
1462 M2 += tmp;
1463 tmp = usr_par.NdNdx1Nw[1]; tmp.Scale(b(0,0)*c(1,0)+b(0,1)*c(1,1));
1464 M2 += tmp;
1465 tmp = usr_par.NdNdx1Nw[2]; tmp.Scale(b(0,0)*c(2,0)+b(0,1)*c(2,1));
1466 M2 += tmp;
1467 tmp = usr_par.NdNdx2Nw[0]; tmp.Scale(b(1,0)*c(0,0)+b(1,1)*c(0,1));
1468 M2 += tmp;
1469 tmp = usr_par.NdNdx2Nw[1]; tmp.Scale(b(1,0)*c(1,0)+b(1,1)*c(1,1));
1470 M2 += tmp;
1471 tmp = usr_par.NdNdx2Nw[2]; tmp.Scale(b(1,0)*c(2,0)+b(1,1)*c(2,1));
1472 M2 += tmp;
1473 M2.Scale(adB/dB);
1474
1475 // Compute the part corresponding to r*u.
1476 tmp = usr_par.NNNw[0]; tmp.Scale(r(0));
1477 M3 += tmp;
1478 tmp = usr_par.NNNw[1]; tmp.Scale(r(1));
1479 M3 += tmp;
1480 tmp = usr_par.NNNw[2]; tmp.Scale(r(2));
1481 M3 += tmp;
1482 M3.Scale(adB);
1483
1484 At += M1;
1485 At += M2;
1486 At += M3;
1487
1488 bt.Multiply('T', 'N', adB, usr_par.NNw, f, 0.0);
1489
1490 return(0);
1491}
1492
1493/* \brief Computes the inverse of a dense matrix.
1494
1495 \param mat [in] - the matrix that is to be inverted
1496 \param inv [in] - the inverse
1497
1498 \return 0 if successful
1499*/
1502{
1503 Epetra_LAPACK lapack;
1504 int dim = mat.M();
1505 int info;
1507 Epetra_SerialDenseVector work(2*dim);
1508
1509 inv.Shape(dim,dim);
1510 inv = mat;
1511
1512 lapack.GETRF(dim, dim, inv.A(), dim, ipiv.A(), &info);
1513 lapack.GETRI(dim, inv.A(), dim, ipiv.A(), work.A(), &dim, &info);
1514
1515 return(0);
1516}
1517
1518/* \brief Computes the determinant of a dense matrix.
1519
1520 \param mat [in] - the matrix
1521
1522 \return the determinant
1523*/
1525{
1526 //Teuchos::LAPACK<int,double> lapack;
1527 Epetra_LAPACK lapack;
1528 double det;
1529 int swaps;
1530 int dim = mat.M();
1531 int info;
1533
1534 Epetra_SerialDenseMatrix mymat(mat);
1535 lapack.GETRF(dim, dim, mymat.A(), dim, ipiv.A(), &info);
1536
1537 det = 1.0;
1538 for (int i=0; i<dim; i++) {
1539 det *= mymat(i,i);
1540 }
1541
1542 swaps = 0;
1543 for (int i=0; i<dim; i++) {
1544 if ((ipiv[i]-1) != i)
1545 swaps++;
1546 }
1547
1548 det *= pow((double)(-1.0),swaps);
1549
1550 return(det);
1551}
1552
1555 Epetra_SerialDenseMatrix & ipcoords,
1557 Epetra_SerialDenseMatrix & pcoords,
1560 const char geomFileBase[],
1561 const bool trace,
1562 const bool dumpAll
1563 )
1564{
1565 int MyPID = Comm.MyPID();
1566
1567 const int FileNameSize = 120;
1568 char FileName[FileNameSize];
1569 TEUCHOS_TEST_FOR_EXCEPT(static_cast<int>(std::strlen(geomFileBase) + 5) > FileNameSize);
1570 sprintf(FileName, "%s.%03d", geomFileBase, MyPID);
1571
1572 {
1573 std::ifstream file_in(FileName);
1574 TEUCHOS_TEST_FOR_EXCEPTION(
1575 file_in.eof(), std::logic_error
1576 ,"Error, the file \""<<FileName<<"\" could not be opened for input!"
1577 );
1578 }
1579
1580 FILE* fid;
1581 fid = fopen(FileName, "r");
1582
1583 if(trace) printf("\nReading node info from %s ...\n", FileName);
1584 int numip = 0, numcp = 0; // # owned nodes, # shared nodes
1585 fscanf(fid, "%d %d", &numip, &numcp);
1586 const int nump = numip + numcp; // # total nodes
1587 if(trace) printf( "\nnumip = %d, numcp = %d, nump = %d\n", numip, numcp, nump );
1588 ipindx.Size(numip);
1589 ipcoords.Shape(numip, 2);
1590 pindx.Size(nump);
1591 pcoords.Shape(nump, 2);
1592 for (int i=0; i<numip; i++) {
1593 fscanf(fid,"%d %lf %lf %*d",&ipindx(i),&ipcoords(i,0),&ipcoords(i,1));
1594 if(trace&&dumpAll) printf("%d %lf %lf\n",ipindx(i),ipcoords(i,0),ipcoords(i,1));
1595 pindx(i) = ipindx(i);
1596 pcoords(i,0) = ipcoords(i,0); pcoords(i,1) = ipcoords(i,1);
1597 }
1598 for (int i=numip; i<nump; i++) {
1599 fscanf(fid,"%d %lf %lf %*d",&pindx(i),&pcoords(i,0),&pcoords(i,1));
1600 }
1601
1602 fscanf(fid,"%*[^\n]"); // Skip to the End of the Line
1603 fscanf(fid,"%*1[\n]"); // Skip One Newline
1604
1605 fscanf(fid,"%*[^\n]"); // Skip to the End of the Line
1606 fscanf(fid,"%*1[\n]"); // Skip One Newline
1607
1608 for (int i=0; i<nump; i++) {
1609 fscanf(fid,"%*[^\n]"); // Skip to the End of the Line
1610 fscanf(fid,"%*1[\n]"); // Skip One Newline
1611 }
1612
1613 if(trace) printf("\nReading element info from %s ...\n", FileName);
1614 int numelems = 0; // # elements on this processor
1615 fscanf(fid, "%d", &numelems);
1616 if(trace) printf( "\nnumelems = %d\n", numelems );
1617 t.Shape(numelems,3);
1618 for (int i=0; i<numelems; i++) {
1619 fscanf(fid, "%d %d %d", &t(i,0), &t(i,1), &t(i,2));
1620 if(trace&&dumpAll) printf("%d %d %d\n", t(i,0), t(i,1), t(i,2));
1621 }
1622
1623 if(trace) printf("\nReading boundry edge info from %s ...\n", FileName);
1624 int numedges = 0; // # boundry edges on this processor
1625 fscanf(fid,"%d",&numedges);
1626 if(trace) printf( "\nnumedges = %d\n", numedges );
1627 e.Shape(numedges,3);
1628 for (int i=0; i<numedges; i++) {
1629 fscanf(fid, "%d %d %d", &e(i,0), &e(i,1), &e(i,2));
1630 if(trace&&dumpAll) printf("%d %d %d\n", e(i,0), e(i,1), e(i,2));
1631 }
1632
1633 fclose(fid);
1634 if(trace) printf("Done reading mesh.\n");
1635
1636 return(0);
1637
1638}
1639
1640/* \brief Performs finite-element assembly of the Hessian of the nonlinear form times an arbitrary vector.
1641
1642 \param Comm [in] - The Epetra (MPI) communicator.
1643 \param ipindx [in] - Vector of NUMIP indices of nodes that are `unique' to a subdomain
1644 (i.e. owned by the corresponding processor).
1645 \param ipcoords [in] - Matrix (NUMIP x 2) of x-y-coordinates of the vertices cooresponding
1646 to indices ipindx: \n
1647 ipcoords(i,0) x-coordinate of node i, \n
1648 ipcoords(i,1) y-coordinate of node i.
1649 \param pindx [in] - Vector of NUMP indices of ALL nodes in a subdomain, including
1650 the shared nodes.
1651 \param pcoords [in] - Matrix (NUMP x 2) of x-y-coordinates of the vertices cooresponding
1652 to indices pindx: \n
1653 pcoords(i,0) x-coordinate of node i, \n
1654 pcoords(i,1) y-coordinate of node i.
1655 \param t [in] - Matrix (ELE x 3) of indices of the vertices in a triangle: \n
1656 t(i,j) index of the j-th vertex in triangle i, where i = 1, ..., ELE
1657 \param y [in] - Reference-counting pointer to the Epetra_MultiVector at which the nonlinear
1658 term is evaluated.
1659 \param s [in] - Reference-counting pointer to the Epetra_MultiVector that is multiplied by the
1660 Hessian of the nonlinear form.
1661 \param lambda [in] - Reference-counting pointer to the Epetra_MultiVector of Lagrange Multipliers.
1662 \param hessvec [out] - Reference-counting pointer to the Epetra_FECrsMatrix containing Hessian of
1663 the nonlinear form times vector product.
1664 \return 0 if successful.
1665
1666 \par Detailed Description:
1667
1668 Assembles the nonlinear term \e hessvec, represented by
1669 \f[
1670 \{N''(y,\lambda)s\}_{j} = \langle g''(y_h) \lambda s,\phi_j \rangle
1671 = \int_{\Omega} g''(y_h(x)) \lambda(x) s(x) \phi_j(x) dx,
1672 \f]
1673 where \f$ g''(y_h) \f$ is given in the local function \e g2pfctn, and \f$\{ \phi_j \}_{j = 1}^{m}\f$ is the
1674 piecewise linear nodal basis for the state space.
1675*/
1677 const Epetra_IntSerialDenseVector & ipindx,
1678 const Epetra_SerialDenseMatrix & ipcoords,
1679 const Epetra_IntSerialDenseVector & pindx,
1680 const Epetra_SerialDenseMatrix & pcoords,
1682 const Teuchos::RCP<const Epetra_MultiVector> & y,
1683 const Teuchos::RCP<const Epetra_MultiVector> & s,
1684 const Teuchos::RCP<const Epetra_MultiVector> & lambda,
1685 Teuchos::RCP<Epetra_FEVector> & hessvec)
1686{
1687
1688 int myPID = Comm.MyPID();
1689 int numProcs = Comm.NumProc();
1690
1691 int numLocNodes = pindx.M();
1692 int numMyLocNodes = ipindx.M();
1693 int numLocElems = t.M();
1694 int numNodesPerElem = 3;
1695
1696 int indexBase = 1;
1697
1698 Epetra_Map standardmap(-1, numMyLocNodes, (int*)ipindx.A(), indexBase, Comm);
1699 Epetra_Map overlapmap(-1, numLocNodes, (int*)pindx.A(), indexBase, Comm);
1700
1701 hessvec = rcp(new Epetra_FEVector(standardmap,1));
1702
1703 int* nodes = new int[numNodesPerElem];
1704 int i=0, j=0, err=0;
1705
1706 // get quadrature nodes and weights
1709 quadrature(2,3,Nodes,Weights);
1710 int numQuadPts = Nodes.M();
1711
1712 // Evaluate nodal basis functions and their derivatives at quadrature points
1713 // N(i,j) = value of the j-th basis function at quadrature node i.
1715 N.Shape(numQuadPts,3);
1716 for (int i=0; i<numQuadPts; i++) {
1717 N(i,0) = 1.0 - Nodes(i,0) - Nodes(i,1);
1718 N(i,1) = Nodes(i,0);
1719 N(i,2) = Nodes(i,1);
1720 }
1721
1722 // Declare quantities needed for the call to the local assembly routine.
1723 Epetra_IntSerialDenseVector epetra_nodes(View, nodes, numNodesPerElem);
1724 Epetra_SerialDenseMatrix vertices(numNodesPerElem, pcoords.N());
1725
1726 Epetra_SerialDenseVector ly; // local entries of y
1727 Epetra_SerialDenseVector Nly; // N*ly
1728 Epetra_SerialDenseVector ls; // local entries of s
1729 Epetra_SerialDenseVector Nls; // N*ls
1730 Epetra_SerialDenseVector llambda; // local entries of lambda
1731 Epetra_SerialDenseVector Nllambda; // N*llambda
1732 Epetra_SerialDenseVector lgfctn; // gfctn(Nly)
1733 Epetra_SerialDenseVector lgfctnNi; // lgfctn.*N(:,i)
1734 Epetra_SerialDenseVector lgfctnNls; // lgfctn.*N(:,i).*llambda.*ls
1735 Epetra_SerialDenseVector lhessvec; // local contribution
1736 // Size and init to zero.
1737 ly.Size(numNodesPerElem);
1738 Nly.Size(numQuadPts);
1739 ls.Size(numNodesPerElem);
1740 Nls.Size(numQuadPts);
1741 llambda.Size(numNodesPerElem);
1742 Nllambda.Size(numQuadPts);
1743 lgfctn.Size(numQuadPts);
1744 lgfctnNi.Size(numQuadPts);
1745 lgfctnNls.Size(numQuadPts);
1746 lhessvec.Size(numNodesPerElem);
1747
1749 double adB;
1750
1751 for(i=0; i<numLocElems; i++) {
1752
1753 nodes[0] = t(i,0); nodes[1] = t(i,1); nodes[2] = t(i,2);
1754 for (j=0; j<numNodesPerElem; j++) {
1755 vertices(j,0) = pcoords(overlapmap.LID(nodes[j]), 0);
1756 vertices(j,1) = pcoords(overlapmap.LID(nodes[j]), 1);
1757 }
1758
1759 // Construct affine transformation matrix.
1760 for(int i=0; i<2; i++) {
1761 B(i,0) = vertices(1,i)-vertices(0,i);
1762 B(i,1) = vertices(2,i)-vertices(0,i);
1763 }
1764 adB = abs(determinant(B));
1765
1766 // Construct local (to each processor) element view of y, s, l.
1767 for (j=0; j<numNodesPerElem; j++) {
1768 ly(j) = (*((*y)(0)))[overlapmap.LID(nodes[j])];
1769 ls(j) = (*((*s)(0)))[overlapmap.LID(nodes[j])];
1770 llambda(j) = (*((*lambda)(0)))[overlapmap.LID(nodes[j])];
1771 }
1772
1773 Nly.Multiply('N', 'N', 1.0, N, ly, 0.0); //N*ly
1774 Nls.Multiply('N', 'N', 1.0, N, ls, 0.0); //N*ls
1775 Nllambda.Multiply('N', 'N', 1.0, N, llambda, 0.0); //N*llambda
1776 g2pfctn(Nly, lgfctn);
1777
1778 for (int i=0; i<numNodesPerElem; i++) {
1779 compproduct(lgfctnNi, lgfctn.A(), N[i]);
1780 compproduct(lgfctnNls, lgfctnNi.A(), Nls.A(), Nllambda.A()); // g2pfctn(Nly).*Nls.*Nllambda.*N(:,i)
1781 lhessvec(i) = adB*lgfctnNls.Dot(Weights);
1782 }
1783
1784 err = hessvec->SumIntoGlobalValues(epetra_nodes, lhessvec);
1785 if (err<0) return(err);
1786 }
1787
1788 // Call global assemble.
1789
1790 err = hessvec->GlobalAssemble();
1791 if (err<0) return(err);
1792
1793 delete [] nodes;
1794
1795 return(0);
1796}
1797
1798/* \brief Componentwise evaluation of the second derivative of the nonlinear reaction term.
1799 \param v [in] - Vector at which the second derivative is evaluated.
1800 \param gv [out] - Vector value.
1801*/
1803 for (int i=0; i<v.M(); i++) {
1804 gv(i) = 6.0*v(i);
1805 }
1806}
1807
1808/* \brief Performs finite-element assembly of the Jacobian of the nonlinear form.
1809
1810 \param Comm [in] - The Epetra (MPI) communicator.
1811 \param ipindx [in] - Vector of NUMIP indices of nodes that are `unique' to a subdomain
1812 (i.e. owned by the corresponding processor).
1813 \param ipcoords [in] - Matrix (NUMIP x 2) of x-y-coordinates of the vertices cooresponding
1814 to indices ipindx: \n
1815 ipcoords(i,0) x-coordinate of node i, \n
1816 ipcoords(i,1) y-coordinate of node i.
1817 \param pindx [in] - Vector of NUMP indices of ALL nodes in a subdomain, including
1818 the shared nodes.
1819 \param pcoords [in] - Matrix (NUMP x 2) of x-y-coordinates of the vertices cooresponding
1820 to indices pindx: \n
1821 pcoords(i,0) x-coordinate of node i, \n
1822 pcoords(i,1) y-coordinate of node i.
1823 \param t [in] - Matrix (ELE x 3) of indices of the vertices in a triangle: \n
1824 t(i,j) index of the j-th vertex in triangle i, where i = 1, ..., ELE
1825 \param y [in] - Reference-counting pointer to the Epetra_MultiVector at which the nonlinear
1826 term is evaluated.
1827 \param Gp [out] - Reference-counting pointer to the Epetra_FECrsMatrix containing the Jacobian
1828 of the nonlinear form.
1829 \return 0 if successful.
1830
1831 \par Detailed Description:
1832
1833 Assembles the nonlinear term \e Gp, represented by
1834 \f[
1835 \{N'(y)\}_{jk} = \langle g'(y_h) \phi_k,\phi_j \rangle = \int_{\Omega} g'(y_h(x)) \phi_k(x) \phi_j(x) dx,
1836 \f]
1837 where \f$ g'(y_h) \f$ is given in the local function \e gpfctn, and \f$\{ \phi_j \}_{j = 1}^{m}\f$ is the
1838 piecewise linear nodal basis for the state space.
1839*/
1841 const Epetra_IntSerialDenseVector & ipindx,
1842 const Epetra_SerialDenseMatrix & ipcoords,
1843 const Epetra_IntSerialDenseVector & pindx,
1844 const Epetra_SerialDenseMatrix & pcoords,
1846 const Teuchos::RCP<const Epetra_MultiVector> & y,
1847 Teuchos::RCP<Epetra_FECrsMatrix> & Gp)
1848{
1849
1850 int myPID = Comm.MyPID();
1851 int numProcs = Comm.NumProc();
1852
1853 int numLocNodes = pindx.M();
1854 int numMyLocNodes = ipindx.M();
1855 int numLocElems = t.M();
1856 int numNodesPerElem = 3;
1857
1858 int indexBase = 1;
1859
1860 Epetra_Map standardmap(-1, numMyLocNodes, (int*)ipindx.A(), indexBase, Comm);
1861 Epetra_Map overlapmap(-1, numLocNodes, (int*)pindx.A(), indexBase, Comm);
1862
1864 Gp = rcp(new Epetra_FECrsMatrix(Copy, standardmap, 0));
1865
1866 int* nodes = new int[numNodesPerElem];
1867 int i=0, j=0, err=0;
1868
1869 // get quadrature nodes and weights
1872 quadrature(2,3,Nodes,Weights);
1873 int numQuadPts = Nodes.M();
1874
1875 // Evaluate nodal basis functions and their derivatives at quadrature points
1876 // N(i,j) = value of the j-th basis function at quadrature node i.
1878 N.Shape(numQuadPts,3);
1879 for (int i=0; i<numQuadPts; i++) {
1880 N(i,0) = 1.0 - Nodes(i,0) - Nodes(i,1);
1881 N(i,1) = Nodes(i,0);
1882 N(i,2) = Nodes(i,1);
1883 }
1884
1885 // Declare quantities needed for the call to the local assembly routine.
1886 Epetra_IntSerialDenseVector epetra_nodes(View, nodes, numNodesPerElem);
1887 Epetra_SerialDenseMatrix vertices(numNodesPerElem, pcoords.N());
1888
1889 Epetra_SerialDenseVector ly; // local entries of y
1890 Epetra_SerialDenseVector Nly; // N*ly
1891 Epetra_SerialDenseVector lgfctn; // gfctn(Nly)
1892 Epetra_SerialDenseVector lgfctnNiNj; // lgfctn.*N(:,i).*N(:,j)
1893 Epetra_SerialDenseMatrix lGp; // local contribution
1894 // Size and init to zero.
1895 ly.Size(numNodesPerElem);
1896 Nly.Size(numQuadPts);
1897 lgfctn.Size(numQuadPts);
1898 lgfctnNiNj.Size(numQuadPts);
1899 lGp.Shape(numNodesPerElem, numNodesPerElem);
1900
1902 double adB;
1903
1904 for(i=0; i<numLocElems; i++) {
1905
1906 nodes[0] = t(i,0); nodes[1] = t(i,1); nodes[2] = t(i,2);
1907 for (j=0; j<numNodesPerElem; j++) {
1908 vertices(j,0) = pcoords(overlapmap.LID(nodes[j]), 0);
1909 vertices(j,1) = pcoords(overlapmap.LID(nodes[j]), 1);
1910 }
1911
1912 // Construct affine transformation matrix.
1913 for(int i=0; i<2; i++) {
1914 B(i,0) = vertices(1,i)-vertices(0,i);
1915 B(i,1) = vertices(2,i)-vertices(0,i);
1916 }
1917 adB = abs(determinant(B));
1918
1919 // Construct local (to each processor) element view of y.
1920 for (j=0; j<numNodesPerElem; j++) {
1921 ly(j) = (*((*y)(0)))[overlapmap.LID(nodes[j])];
1922 }
1923
1924 Nly.Multiply('N', 'N', 1.0, N, ly, 0.0);
1925 gpfctn(Nly, lgfctn);
1926
1927 for (int i=0; i<numNodesPerElem; i++) {
1928 for (int j=0; j<numNodesPerElem; j++) {
1929 compproduct(lgfctnNiNj, lgfctn.A(), N[i], N[j]);
1930 lGp(i,j) = adB*lgfctnNiNj.Dot(Weights);
1931 }
1932 }
1933
1934 err = Gp->InsertGlobalValues(epetra_nodes, lGp, format);
1935 if (err<0) return(err);
1936 }
1937
1938 // Call global assemble.
1939
1940 err = Gp->GlobalAssemble();
1941 if (err<0) return(err);
1942
1943 delete [] nodes;
1944
1945 return(0);
1946}
1947
1948/* \brief Componentwise evaluation of the first derivative of the nonlinear reaction term.
1949 \param v [in] - Vector at which the first derivative is evaluated.
1950 \param gv [out] - Vector value.
1951*/
1953 for (int i=0; i<v.M(); i++) {
1954 gv(i) = 3.0*pow(v(i),2)-1.0;
1955 }
1956}
1957
1958/* \brief Performs finite-element assembly of the nonlinear reaction term.
1959
1960 \param Comm [in] - The Epetra (MPI) communicator.
1961 \param ipindx [in] - Vector of NUMIP indices of nodes that are `unique' to a subdomain
1962 (i.e. owned by the corresponding processor).
1963 \param ipcoords [in] - Matrix (NUMIP x 2) of x-y-coordinates of the vertices cooresponding
1964 to indices ipindx: \n
1965 ipcoords(i,0) x-coordinate of node i, \n
1966 ipcoords(i,1) y-coordinate of node i.
1967 \param pindx [in] - Vector of NUMP indices of ALL nodes in a subdomain, including
1968 the shared nodes.
1969 \param pcoords [in] - Matrix (NUMP x 2) of x-y-coordinates of the vertices cooresponding
1970 to indices pindx: \n
1971 pcoords(i,0) x-coordinate of node i, \n
1972 pcoords(i,1) y-coordinate of node i.
1973 \param t [in] - Matrix (ELE x 3) of indices of the vertices in a triangle: \n
1974 t(i,j) index of the j-th vertex in triangle i, where i = 1, ..., ELE
1975 \param y [out] - Reference-counting pointer to the Epetra_MultiVector at which the nonlinear
1976 form is evaluated.
1977 \param g [out] - Reference-counting pointer to the Epetra_FEVector containing the value
1978 of the nonlinear form.
1979 \return 0 if successful.
1980
1981 \par Detailed Description:
1982
1983 Assembles the nonlinear term \e g, represented by
1984 \f[
1985 \{N(y)\}_{j} = \langle g(y_h),\phi_j \rangle = \int_{\Omega} g(y_h(x)) \phi_j(x) dx,
1986 \f]
1987 where \f$ g(y_h) \f$ is given in the local function \e gfctn, and \f$\{ \phi_j \}_{j = 1}^{m}\f$ is the
1988 piecewise linear nodal basis for the state space.
1989*/
1991 const Epetra_IntSerialDenseVector & ipindx,
1992 const Epetra_SerialDenseMatrix & ipcoords,
1993 const Epetra_IntSerialDenseVector & pindx,
1994 const Epetra_SerialDenseMatrix & pcoords,
1996 const Teuchos::RCP<const Epetra_MultiVector> & y,
1997 Teuchos::RCP<Epetra_FEVector> & g)
1998{
1999
2000 int myPID = Comm.MyPID();
2001 int numProcs = Comm.NumProc();
2002
2003 int numLocNodes = pindx.M();
2004 int numMyLocNodes = ipindx.M();
2005 int numLocElems = t.M();
2006 int numNodesPerElem = 3;
2007
2008 int indexBase = 1;
2009
2010 Epetra_Map standardmap(-1, numMyLocNodes, (int*)ipindx.A(), indexBase, Comm);
2011 Epetra_Map overlapmap(-1, numLocNodes, (int*)pindx.A(), indexBase, Comm);
2012
2013 g = rcp(new Epetra_FEVector(standardmap,1));
2014
2015 int* nodes = new int[numNodesPerElem];
2016 int i=0, j=0, err=0;
2017
2018 // get quadrature nodes and weights
2021 quadrature(2,3,Nodes,Weights);
2022 int numQuadPts = Nodes.M();
2023
2024 // Evaluate nodal basis functions and their derivatives at quadrature points
2025 // N(i,j) = value of the j-th basis function at quadrature node i.
2027 N.Shape(numQuadPts,3);
2028 for (int i=0; i<numQuadPts; i++) {
2029 N(i,0) = 1.0 - Nodes(i,0) - Nodes(i,1);
2030 N(i,1) = Nodes(i,0);
2031 N(i,2) = Nodes(i,1);
2032 }
2033
2034 // Declare quantities needed for the call to the local assembly routine.
2035 Epetra_IntSerialDenseVector epetra_nodes(View, nodes, numNodesPerElem);
2036 Epetra_SerialDenseMatrix vertices(numNodesPerElem, pcoords.N());
2037
2038 Epetra_SerialDenseVector ly; // local entries of y
2039 Epetra_SerialDenseVector Nly; // N*ly
2040 Epetra_SerialDenseVector lgfctn; // gfctn(Nly)
2041 Epetra_SerialDenseVector lgfctnNi; // lgfctn.*N(:,i)
2042 Epetra_SerialDenseVector lg; // local contribution
2043 // Size and init to zero.
2044 ly.Size(numNodesPerElem);
2045 Nly.Size(numQuadPts);
2046 lgfctn.Size(numQuadPts);
2047 lgfctnNi.Size(numQuadPts);
2048 lg.Size(numNodesPerElem);
2049
2051 double adB;
2052
2053 for(i=0; i<numLocElems; i++) {
2054
2055 nodes[0] = t(i,0); nodes[1] = t(i,1); nodes[2] = t(i,2);
2056 for (j=0; j<numNodesPerElem; j++) {
2057 vertices(j,0) = pcoords(overlapmap.LID(nodes[j]), 0);
2058 vertices(j,1) = pcoords(overlapmap.LID(nodes[j]), 1);
2059 }
2060
2061 // Construct affine transformation matrix.
2062 for(int i=0; i<2; i++) {
2063 B(i,0) = vertices(1,i)-vertices(0,i);
2064 B(i,1) = vertices(2,i)-vertices(0,i);
2065 }
2066 adB = abs(determinant(B));
2067
2068 // Construct local (to each processor) element view of y.
2069 for (j=0; j<numNodesPerElem; j++) {
2070 ly(j) = (*((*y)(0)))[overlapmap.LID(nodes[j])];
2071 }
2072
2073 Nly.Multiply('N', 'N', 1.0, N, ly, 0.0);
2074 gfctn(Nly, lgfctn);
2075
2076 for (int i=0; i<numNodesPerElem; i++) {
2077 compproduct(lgfctnNi, lgfctn.A(), N[i]);
2078 lg(i) = adB*lgfctnNi.Dot(Weights);
2079 }
2080
2081 err = g->SumIntoGlobalValues(epetra_nodes, lg);
2082 if (err<0) return(err);
2083 }
2084
2085 // Call global assemble.
2086
2087 err = g->GlobalAssemble();
2088 if (err<0) return(err);
2089
2090 delete [] nodes;
2091
2092 return(0);
2093}
2094
2095
2096/* \brief Componentwise evaluation of the nonlinear reaction term.
2097 \param v [in] - Vector at which the nonlinear function is evaluated.
2098 \param gv [out] - Vector value.
2099*/
2101 for (int i=0; i<v.M(); i++) {
2102 gv(i) = pow(v(i),3)-v(i);
2103 }
2104}
2105
2106/* ======== ================ *
2107 * function CrsMatrix2MATLAB *
2108 * ======== ================ *
2109 *
2110 * Print out a CrsMatrix in a MATLAB format. Each processor prints out
2111 * its part, starting from proc 0 to proc NumProc-1. The first line of
2112 * each processor's output states the number of local rows and of
2113 * local nonzero elements.
2114 *
2115 *
2116 * Return code: true if matrix has been printed out
2117 * ----------- false otherwise
2118 *
2119 * Parameters:
2120 * ----------
2121 *
2122 * - Epetra_CrsMatrix reference to the distributed CrsMatrix to
2123 * print out
2124 * - std::ostream & reference to output stream
2125 */
2126
2127bool GLpApp::CrsMatrix2MATLAB(const Epetra_CrsMatrix & A, std::ostream & outfile)
2128{
2129
2130 int MyPID = A.Comm().MyPID();
2131 int NumProc = A.Comm().NumProc();
2132
2133 // work only on transformed matrices;
2134 if( A.IndicesAreLocal() == false ) {
2135 if( MyPID == 0 ) {
2136 std::cerr << "ERROR in "<< __FILE__ << ", line " << __LINE__ << std::endl;
2137 std::cerr << "Function CrsMatrix2MATLAB accepts\n";
2138 std::cerr << "transformed matrices ONLY. Please call A.TransformToLoca()\n";
2139 std::cerr << "on your matrix A to that purpose.\n";
2140 std::cerr << "Now returning...\n";
2141 }
2142 return false;
2143 }
2144
2145 int NumMyRows = A.NumMyRows(); // number of rows on this process
2146 int NumNzRow; // number of nonzero elements for each row
2147 int NumEntries; // number of extracted elements for each row
2148 int NumGlobalRows; // global dimensio of the problem
2149 int GlobalRow; // row in global ordering
2150 int NumGlobalNonzeros; // global number of nonzero elements
2151
2152 NumGlobalRows = A.NumGlobalRows();
2153 NumGlobalNonzeros = A.NumGlobalNonzeros();
2154
2155 // print out on std::cout if no filename is provided
2156
2157 int IndexBase = A.IndexBase(); // MATLAB starts from 0
2158 if( IndexBase == 0 )
2159 IndexBase = 1;
2160 else if ( IndexBase == 1)
2161 IndexBase = 0;
2162
2163 // write on file the dimension of the matrix
2164
2165 if( MyPID==0 ) {
2166 outfile << "A = spalloc(";
2167 outfile << NumGlobalRows << ',' << NumGlobalRows;
2168 outfile << ',' << NumGlobalNonzeros << ");\n";
2169 }
2170
2171 for( int Proc=0 ; Proc<NumProc ; ++Proc ) {
2172 A.Comm().Barrier();
2173 if( MyPID == Proc ) {
2174
2175 outfile << "\n\n% On proc " << Proc << ": ";
2176 outfile << NumMyRows << " rows and ";
2177 outfile << A.NumMyNonzeros() << " nonzeros\n";
2178
2179 // cycle over all local rows to find out nonzero elements
2180 for( int MyRow=0 ; MyRow<NumMyRows ; ++MyRow ) {
2181
2182 GlobalRow = A.GRID(MyRow);
2183
2184 NumNzRow = A.NumMyEntries(MyRow);
2185 double *Values = new double[NumNzRow];
2186 int *Indices = new int[NumNzRow];
2187
2188 A.ExtractMyRowCopy(MyRow, NumNzRow,
2189 NumEntries, Values, Indices);
2190 // print out the elements with MATLAB syntax
2191 for( int j=0 ; j<NumEntries ; ++j ) {
2192 outfile << "A(" << GlobalRow + IndexBase
2193 << "," << A.GCID(Indices[j]) + IndexBase
2194 << ") = " << Values[j] << ";\n";
2195 }
2196
2197 delete Values;
2198 delete Indices;
2199 }
2200
2201 }
2202 A.Comm().Barrier();
2203 }
2204
2205 return true;
2206
2207}
2208
2209
2210/* ======== ============= *
2211 * function Vector2MATLAB *
2212 * ======== ============= *
2213 *
2214 * Print out a Epetra_Vector in a MATLAB format. Each processor prints out
2215 * its part, starting from proc 0 to proc NumProc-1. The first line of
2216 * each processor's output states the number of local rows and of
2217 * local nonzero elements.
2218 *
2219 * Return code: true if vector has been printed out
2220 * ----------- false otherwise
2221 *
2222 * Parameters:
2223 * ----------
2224 *
2225 * - Epetra_Vector reference to vector
2226 * - std::ostream & reference to output stream
2227 */
2228
2229bool GLpApp::Vector2MATLAB( const Epetra_Vector & v, std::ostream & outfile)
2230{
2231
2232 int MyPID = v.Comm().MyPID();
2233 int NumProc = v.Comm().NumProc();
2234 int MyLength = v.MyLength();
2235 int GlobalLength = v.GlobalLength();
2236
2237 // print out on std::cout if no filename is provided
2238
2239 // write on file the dimension of the matrix
2240
2241 if( MyPID == 0 ) outfile << "v = zeros(" << GlobalLength << ",1)\n";
2242
2243 int NumMyElements = v.Map().NumMyElements();
2244 // get update list
2245 int * MyGlobalElements = v.Map().MyGlobalElements( );
2246
2247 int Row;
2248
2249 int IndexBase = v.Map().IndexBase(); // MATLAB starts from 0
2250 if( IndexBase == 0 )
2251 IndexBase = 1;
2252 else if ( IndexBase == 1)
2253 IndexBase = 0;
2254
2255 for( int Proc=0 ; Proc<NumProc ; ++Proc ) {
2256 v.Comm().Barrier();
2257 if( MyPID == Proc ) {
2258
2259 outfile << "% On proc " << Proc << ": ";
2260 outfile << MyLength << " rows of ";
2261 outfile << GlobalLength << " elements\n";
2262
2263 for( Row=0 ; Row<MyLength ; ++Row ) {
2264 outfile << "v(" << MyGlobalElements[Row] + IndexBase
2265 << ") = " << v[Row] << ";\n";
2266 }
2267
2268 }
2269
2270 v.Comm().Barrier();
2271 }
2272
2273 return true;
2274
2275} /* Vector2MATLAB */
2276
2277
2278/* ======== =============== *
2279 * function FEVector2MATLAB *
2280 * ======== =============== *
2281 *
2282 * Print out a Epetra_Vector in a MATLAB format. Each processor prints out
2283 * its part, starting from proc 0 to proc NumProc-1. The first line of
2284 * each processor's output states the number of local rows and of
2285 * local nonzero elements.
2286 *
2287 * Return code: true if vector has been printed out
2288 * ----------- false otherwise
2289 *
2290 * Parameters:
2291 * ----------
2292 *
2293 * - Epetra_FEVector reference to FE vector
2294 * - std::ostream & reference to output stream
2295 */
2296
2297bool GLpApp::FEVector2MATLAB( const Epetra_FEVector & v, std::ostream & outfile)
2298{
2299
2300 int MyPID = v.Comm().MyPID();
2301 int NumProc = v.Comm().NumProc();
2302 int MyLength = v.MyLength();
2303 int GlobalLength = v.GlobalLength();
2304
2305 // print out on std::cout if no filename is provided
2306
2307 // write on file the dimension of the matrix
2308
2309 if( MyPID == 0 ) outfile << "v = zeros(" << GlobalLength << ",1)\n";
2310
2311 int NumMyElements = v.Map().NumMyElements();
2312 // get update list
2313 int * MyGlobalElements = v.Map().MyGlobalElements( );
2314
2315 int Row;
2316
2317 int IndexBase = v.Map().IndexBase(); // MATLAB starts from 0
2318 if( IndexBase == 0 )
2319 IndexBase = 1;
2320 else if ( IndexBase == 1)
2321 IndexBase = 0;
2322
2323 for( int Proc=0 ; Proc<NumProc ; ++Proc ) {
2324 v.Comm().Barrier();
2325 if( MyPID == Proc ) {
2326
2327 outfile << "% On proc " << Proc << ": ";
2328 outfile << MyLength << " rows of ";
2329 outfile << GlobalLength << " elements\n";
2330
2331 for( Row=0 ; Row<MyLength ; ++Row ) {
2332 outfile << "v(" << MyGlobalElements[Row] + IndexBase
2333 << ") = " << v[0][Row] << ";\n";
2334 }
2335
2336 }
2337
2338 v.Comm().Barrier();
2339 }
2340
2341 return true;
2342
2343} /* FEVector2MATLAB */
2344
2345
2346/* \brief Returns the nodes and weights for the integration \n
2347 on the interval [0,1] (dim = 1) \n
2348 on the triangle with vertices (0,0), (1,0), (0,1) (if dim = 2) \n
2349 on the tetrahedron with vertices (0,0,0), (1,0,0), (0,1,0), (0,0,1) (if dim = 3).
2350
2351 \param dim [in] - spatial dimension (dim = 1, 2)
2352 \param order [in] - required degree of polynomials that integrate exactly
2353 \param nodes [out] - Matrix in which the i-th row of nodes gives coordinates of the i-th quadrature node
2354 \param weights [out] - quadrature weights
2355
2356 \return 0 if successful
2357*/
2358int GLpApp::quadrature(const int dim, const int order,
2360 Epetra_SerialDenseVector & weights)
2361{
2362
2363 if (dim == 1) {
2364
2365 // Gauss quadrature nodes and weights on the interval [0,1]
2366
2367 if (order == 1) {
2368 nodes.Shape(1,1);
2369 nodes(0,0) = 0.5;
2370 weights.Size(1);
2371 weights(0) = 1.0;
2372 }
2373 else if (order == 2) {
2374 nodes.Shape(2,1);
2375 nodes(0,0) = (1.0-1.0/sqrt(3.0))/2.0;
2376 nodes(1,0) = (1.0+1.0/sqrt(3.0))/2.0;
2377 weights.Size(2);
2378 weights(0) = 0.5;
2379 weights(1) = 0.5;
2380 }
2381 else if (order == 3) {
2382 nodes.Shape(3,1);
2383 nodes(0,0) = (1.0-sqrt(3.0/5.0))/2.0;
2384 nodes(1,0) = 0.5;
2385 nodes(2,0) = (1.0+sqrt(3.0/5.0))/2.0;
2386 weights.Size(3);
2387 weights(0) = 5.0/18.0;
2388 weights(1) = 4.0/9.0;
2389 weights(2) = 5.0/18.0;
2390 }
2391 else {
2392 std::cout << "Quadrature for dim = " << dim << " and order = ";
2393 std::cout << order << " not available.\n";
2394 exit(-1);
2395 }
2396
2397 }
2398 else if (dim == 2) {
2399
2400 // Quadrature nodes and weights on the unit simplex with
2401 // vertices (0,0), (1,0), and (0,1).
2402
2403 if (order == 1) {
2404 nodes.Shape(1,2);
2405 nodes(0,0) = 1.0/3.0; nodes (0,1) = 1.0/3.0;
2406 weights.Size(1);
2407 weights(0) = 0.5;
2408 }
2409 else if (order == 2) {
2410 nodes.Shape(3,2);
2411 nodes(0,0) = 1.0/6.0; nodes (0,1) = 1.0/6.0;
2412 nodes(1,0) = 2.0/3.0; nodes (1,1) = 1.0/6.0;
2413 nodes(2,0) = 1.0/6.0; nodes (2,1) = 2.0/3.0;
2414 weights.Size(3);
2415 weights(0) = 1.0/6.0;
2416 weights(1) = 1.0/6.0;
2417 weights(2) = 1.0/6.0;
2418 }
2419 else if (order == 3) {
2420 nodes.Shape(4,2);
2421 nodes(0,0) = 1.0/3.0; nodes (0,1) = 1.0/3.0;
2422 nodes(1,0) = 3.0/5.0; nodes (1,1) = 1.0/5.0;
2423 nodes(2,0) = 1.0/5.0; nodes (2,1) = 3.0/5.0;
2424 nodes(3,0) = 1.0/5.0; nodes (3,1) = 1.0/5.0;
2425 weights.Size(4);
2426 weights(0) = -9.0/32.0;
2427 weights(1) = 25.0/96.0;
2428 weights(2) = 25.0/96.0;
2429 weights(3) = 25.0/96.0;
2430 }
2431 else {
2432 std::cout << "Quadrature for dim = " << dim << " and order = ";
2433 std::cout << order << " not available.\n";
2434 exit(-1);
2435 }
2436
2437 }
2438 else {
2439 std::cout << "Quadrature for dim = " << dim << " not available.\n";
2440 exit(-1);
2441 }
2442
2443 return(0);
2444}
Insert
View
Copy
Transform to form the explicit transpose of a Epetra_RowMatrix.
float DOT(const int N, const float *X, const float *Y, const int INCX=1, const int INCY=1) const
float NRM2(const int N, const float *X, const int INCX=1) const
int MyGlobalElements(int *MyGlobalElementList) const
virtual void Print(std::ostream &os) const
int LID(int GID) const
int IndexBase() const
int NumGlobalElements() const
int NumMyElements() const
virtual int NumProc() const=0
virtual int MyPID() const=0
virtual void Barrier() const=0
int GCID(int LCID_in) const
int NumGlobalNonzeros() const
int NumMyNonzeros() const
int ExtractGlobalRowCopy(int_type Row, int Length, int &NumEntries, double *values, int_type *Indices) const
int GRID(int LRID_in) const
const Epetra_Comm & Comm() const
int NumGlobalRows() const
int IndexBase() const
bool IndicesAreLocal() const
int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
int NumMyEntries(int Row) const
int NumMyRows() const
const Epetra_BlockMap & Map() const
const Epetra_Comm & Comm() const
int Shape(int NumRows, int NumCols)
const int * A() const
int Resize(int Length_in)
int Size(int Length_in)
void GETRF(const int M, const int N, float *A, const int LDA, int *IPIV, int *INFO) const
void GETRI(const int N, float *A, const int LDA, int *IPIV, float *WORK, const int *LWORK, int *INFO) const
int GlobalLength() const
int MyLength() const
double * A() const
virtual void Print(std::ostream &os) const
int Scale(double ScalarA)
int Shape(int NumRows, int NumCols)
int Multiply(char TransA, char TransB, double ScalarAB, const Epetra_SerialDenseMatrix &A, const Epetra_SerialDenseMatrix &B, double ScalarThis)
virtual void Print(std::ostream &os) const
int Size(int Length_in)
int Resize(int Length_in)
double * Values() const
double Dot(const Epetra_SerialDenseVector &x) const
Teuchos::RCP< Epetra_FEVector > getNy()
Teuchos::RCP< Epetra_FECrsMatrix > H_
Volume mass matrix.
Teuchos::RCP< Epetra_FECrsMatrix > R_
Edge mass matrix.
Teuchos::RCP< Epetra_FECrsMatrix > getB()
Teuchos::RCP< Epetra_FECrsMatrix > getR()
Teuchos::RCP< Epetra_IntSerialDenseVector > ipindx_
Global nodes (interior, nonoverlapping) in this subdomain.
Teuchos::RCP< const Epetra_SerialDenseMatrix > getpcoords()
Teuchos::RCP< Epetra_SerialDenseMatrix > pcoords_
Coordinates of all nodes in this subdomain.
Teuchos::RCP< Epetra_FEVector > q_
The desired state.
Teuchos::RCP< const Epetra_SerialDenseMatrix > getipcoords()
Teuchos::RCP< const Epetra_IntSerialDenseVector > getipindx()
Teuchos::RCP< Epetra_FEVector > b_
Right-hand side of the PDE.
int solveAugsys(const Teuchos::RCP< const Epetra_MultiVector > &rhsy, const Teuchos::RCP< const Epetra_MultiVector > &rhsu, const Teuchos::RCP< const Epetra_MultiVector > &rhsp, const Teuchos::RCP< Epetra_MultiVector > &y, const Teuchos::RCP< Epetra_MultiVector > &u, const Teuchos::RCP< Epetra_MultiVector > &p, double tol)
Solves augmented system.
GLpYUEpetraDataPool(Teuchos::RCP< const Epetra_Comm > const &commptr, const double beta, const double len_x, const double len_y, const int local_nx, const int local_ny, const char myfile[], const bool trace)
Teuchos::RCP< const Epetra_Comm > getCommPtr()
Teuchos::RCP< Epetra_FEVector > Ny_
Teuchos::RCP< Epetra_FEVector > getb()
void computeNy(const Teuchos::RCP< const Epetra_MultiVector > &y)
Calls the function that computes the nonlinear term.
Teuchos::RCP< Epetra_CrsMatrix > Augmat_
Augmented system matrix: [ I Jac* ] [Jac 0 ].
void computeNpy(const Teuchos::RCP< const Epetra_MultiVector > &y)
Calls the function that computes the Jacobian of the nonlinear term.
Teuchos::RCP< const Epetra_IntSerialDenseMatrix > gete()
void computeAll(const GenSQP::Vector &x)
Calls functions to compute nonlinear quantities and the augmented system matrix.
Teuchos::RCP< Epetra_IntSerialDenseMatrix > t_
Elements (this includes all overlapping nodes).
Teuchos::RCP< Epetra_FECrsMatrix > getA()
void PrintVec(const Teuchos::RCP< const Epetra_Vector > &x)
Outputs the solution vector to files.
Teuchos::RCP< Epetra_CrsMatrix > getAugmat()
Teuchos::RCP< Epetra_FECrsMatrix > A_
Volume stiffness matrix.
void computeAugmat()
Assembles the augmented system (KKT-type) matrix.
Teuchos::RCP< Epetra_SerialDenseMatrix > ipcoords_
Coordinates of nodes that are unique to this subdomain.
Teuchos::RCP< Epetra_FECrsMatrix > Npy_
Jacobian of the nonlinear term.
Teuchos::RCP< Epetra_IntSerialDenseMatrix > e_
Edges.
Teuchos::RCP< Epetra_FECrsMatrix > getNpy()
Teuchos::RCP< const Epetra_IntSerialDenseMatrix > gett()
Teuchos::RCP< Epetra_FEVector > getq()
Teuchos::RCP< Epetra_IntSerialDenseVector > pindx_
Global nodes (interior + shared, overlapping) in this subdomain.
Teuchos::RCP< Epetra_FECrsMatrix > B_
Control/state mass matrix.
double beta_
Regularization parameter.
Teuchos::RCP< const Epetra_IntSerialDenseVector > getpindx()
Teuchos::RCP< const Epetra_Comm > commptr_
Teuchos::RCP< Epetra_FECrsMatrix > getH()
void Print(std::ostream &os) const
Epetra_SerialDenseMatrix S3
Epetra_SerialDenseMatrix Nodes
Epetra_SerialDenseVector Nw
Epetra_SerialDenseMatrix S2
Epetra_SerialDenseMatrix * NdNdx2Nw
Epetra_SerialDenseMatrix Nx1
Epetra_SerialDenseMatrix * NdNdx1Nw
Epetra_SerialDenseMatrix N
Epetra_SerialDenseMatrix S1
Epetra_SerialDenseMatrix Nx2
Epetra_SerialDenseMatrix * NNNw
Epetra_SerialDenseVector Weights
Epetra_SerialDenseMatrix NNw
Provides the interface to generic abstract vector libraries.
const int N
bool CrsMatrix2MATLAB(const Epetra_CrsMatrix &, std::ostream &)
void gpfctn(const Epetra_SerialDenseVector &v, Epetra_SerialDenseVector &gv)
bool FEVector2MATLAB(const Epetra_FEVector &, std::ostream &)
int quadrature(const int, const int, Epetra_SerialDenseMatrix &, Epetra_SerialDenseVector &)
void rect2DMeshGenerator(const int numProc, const int procRank, const double len_x, const double len_y, const int local_nx, const int local_ny, const int bndy_marker, Epetra_IntSerialDenseVector *ipindx_out, Epetra_SerialDenseMatrix *ipcoords_out, Epetra_IntSerialDenseVector *pindx_out, Epetra_SerialDenseMatrix *pcoords_out, Epetra_IntSerialDenseMatrix *t_out, Epetra_IntSerialDenseMatrix *e_out, std::ostream *out, const bool dumpAll)
Generate a simple rectangular 2D triangular mesh that is only partitioned between processors in the y...
int assemblyFECrs(const Epetra_Comm &, const Epetra_IntSerialDenseVector &, const Epetra_SerialDenseMatrix &, const Epetra_IntSerialDenseVector &, const Epetra_SerialDenseMatrix &, const Epetra_IntSerialDenseMatrix &, const Epetra_IntSerialDenseMatrix &, Teuchos::RCP< Epetra_FECrsMatrix > &, Teuchos::RCP< Epetra_FEVector > &)
int meshreader(const Epetra_Comm &, Epetra_IntSerialDenseVector &, Epetra_SerialDenseMatrix &, Epetra_IntSerialDenseVector &, Epetra_SerialDenseMatrix &, Epetra_IntSerialDenseMatrix &, Epetra_IntSerialDenseMatrix &, const char geomFileBase[], const bool trace=true, const bool dumpAll=false)
int assemble(const Epetra_Comm &, const Epetra_IntSerialDenseVector &, const Epetra_SerialDenseMatrix &, const Epetra_IntSerialDenseVector &, const Epetra_SerialDenseMatrix &, const Epetra_IntSerialDenseMatrix &, const Epetra_IntSerialDenseMatrix &, Teuchos::RCP< Epetra_FECrsMatrix > &, Teuchos::RCP< Epetra_FECrsMatrix > &, Teuchos::RCP< Epetra_FEVector > &)
int inverse(const Epetra_SerialDenseMatrix &, Epetra_SerialDenseMatrix &)
double determinant(const Epetra_SerialDenseMatrix &)
int nonlinjac(const Epetra_Comm &, const Epetra_IntSerialDenseVector &, const Epetra_SerialDenseMatrix &, const Epetra_IntSerialDenseVector &, const Epetra_SerialDenseMatrix &, const Epetra_IntSerialDenseMatrix &, const Teuchos::RCP< const Epetra_MultiVector > &, Teuchos::RCP< Epetra_FECrsMatrix > &)
int lassembly(const Epetra_SerialDenseMatrix &, const Epetra_SerialDenseVector &, const Epetra_SerialDenseMatrix &, const Epetra_SerialDenseVector &, const Epetra_SerialDenseVector &, const Usr_Par &, Epetra_SerialDenseMatrix &, Epetra_SerialDenseVector &)
int nonlinvec(const Epetra_Comm &, const Epetra_IntSerialDenseVector &, const Epetra_SerialDenseMatrix &, const Epetra_IntSerialDenseVector &, const Epetra_SerialDenseMatrix &, const Epetra_IntSerialDenseMatrix &, const Teuchos::RCP< const Epetra_MultiVector > &, Teuchos::RCP< Epetra_FEVector > &)
int assemble_bdry(const Epetra_Comm &Comm, const Epetra_IntSerialDenseVector &ipindx, const Epetra_SerialDenseMatrix &ipcoords, const Epetra_IntSerialDenseVector &pindx, const Epetra_SerialDenseMatrix &pcoords, const Epetra_IntSerialDenseMatrix &t, const Epetra_IntSerialDenseMatrix &e, Teuchos::RCP< Epetra_FECrsMatrix > *B, Teuchos::RCP< Epetra_FECrsMatrix > *R)
bool Vector2MATLAB(const Epetra_Vector &, std::ostream &)
int compproduct(Epetra_SerialDenseVector &, double *, double *)
void g2pfctn(const Epetra_SerialDenseVector &, Epetra_SerialDenseVector &)
void gfctn(const Epetra_SerialDenseVector &, Epetra_SerialDenseVector &)
int nonlinhessvec(const Epetra_Comm &, const Epetra_IntSerialDenseVector &, const Epetra_SerialDenseMatrix &, const Epetra_IntSerialDenseVector &, const Epetra_SerialDenseMatrix &, const Epetra_IntSerialDenseMatrix &, const Teuchos::RCP< const Epetra_MultiVector > &, const Teuchos::RCP< const Epetra_MultiVector > &, const Teuchos::RCP< const Epetra_MultiVector > &, Teuchos::RCP< Epetra_FEVector > &)
std::ostream & operator<<(std::ostream &, const Usr_Par &)