EpetraExt Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
50 #include "rect2DMeshGenerator.hpp"
51 
52 #include "Amesos_Klu.h"
53 #include "EpetraExt_MatrixMatrix.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"
69 #include "Teuchos_RCP.hpp"
70 #include "Teuchos_Assert.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)
82 using namespace std;
83 
84 
85 // Define to see all debug output for mesh generation
86 //#define GLPYUEPETRA_DATAPOOL_DUMP_ALL_MESH
87 
88 
89 namespace GLpApp {
90 
91 //
92 // Declarations
93 //
94 
95 const double GLp_pi = 3.14159265358979323846;
96 
97 std::ostream& operator <<(std::ostream &, const Usr_Par &);
98 
99 bool CrsMatrix2MATLAB(const Epetra_CrsMatrix &, std::ostream &);
100 
101 bool Vector2MATLAB( const Epetra_Vector &, std::ostream &);
102 
103 bool FEVector2MATLAB( const Epetra_FEVector &, std::ostream &);
104 
105 int quadrature(const int, const int, Epetra_SerialDenseMatrix &,
107 
108 int meshreader(
109  const Epetra_Comm &,
116  const char geomFileBase[],
117  const bool trace = true,
118  const bool dumpAll = false
119  );
120 
122  const Epetra_SerialDenseVector &,
123  const Epetra_SerialDenseMatrix &,
124  const Epetra_SerialDenseVector &,
125  const Epetra_SerialDenseVector &,
126  const Usr_Par &,
129 
130 int assemblyFECrs(const Epetra_Comm &,
132  const Epetra_SerialDenseMatrix &,
134  const Epetra_SerialDenseMatrix &,
139 
140 int assemblyFECrs(const Epetra_Comm &,
142  const Epetra_SerialDenseMatrix &,
144  const Epetra_SerialDenseMatrix &,
149  bool);
150 
151 int assemble(const Epetra_Comm &,
153  const Epetra_SerialDenseMatrix &,
155  const Epetra_SerialDenseMatrix &,
161 
162 int 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
172  );
173 
174 int nonlinvec(const Epetra_Comm &,
176  const Epetra_SerialDenseMatrix &,
178  const Epetra_SerialDenseMatrix &,
182 
183 int nonlinjac(const Epetra_Comm &,
185  const Epetra_SerialDenseMatrix &,
187  const Epetra_SerialDenseMatrix &,
191 
192 int nonlinhessvec(const Epetra_Comm &,
194  const Epetra_SerialDenseMatrix &,
196  const Epetra_SerialDenseMatrix &,
202 
203 int compproduct(Epetra_SerialDenseVector &, double *, double *);
204 
205 int compproduct(Epetra_SerialDenseVector &, double *, double *, double *);
206 
207 double determinant(const Epetra_SerialDenseMatrix &);
208 
210 
211 int quadrature(
212  const int, const int, Epetra_SerialDenseMatrix &,
214 
216 
218 
220 
221 //
222 // GLpYUEpetraDataPool
223 //
224 
225 GLpYUEpetraDataPool::GLpYUEpetraDataPool(
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 {
244 
245  if( myfile && myfile[0] != '\0' ) {
246  meshreader(*commptr_,*ipindx_,*ipcoords_,*pindx_,*pcoords_,*t_,*e_,myfile,trace);
247  }
248  else {
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.
282  (Teuchos::dyn_cast<GenSQP::YUEpetraVector>(const_cast<GenSQP::Vector&>(x))).getYVector();
283 
284  computeNy(ey);
285 
286  computeNpy(ey);
287 
288  computeAugmat();
289 
290 }
291 
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 */
492  return 0;
493 }
494 
496 
498 
500 
502 
504 
506 
508 
510 
512 
514 
516 
518 
520 
522 
524 
526 
528 
529 
531 {
532  Epetra_Map overlapmap(-1, pindx_->M(), (int*)(pindx_)->A(), 1, *commptr_);
533  Epetra_Map standardmap(A_->DomainMap());
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 
542 {
543  Epetra_Map overlapmap(-1, pindx_->M(), (int*)(pindx_)->A(), 1, *commptr_);
544  Epetra_Map standardmap(A_->DomainMap());
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 
683 {
684  Vector2MATLAB(*x, std::cout);
685 }
686 
688  Epetra_BLAS blas;
689  Epetra_SerialDenseVector product(4);
690 
691  // get nodes and weights
692  quadrature(2,3,Nodes,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 
783 void 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 
816 std::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
900  )
901 {
902 
903  using Teuchos::rcp;
904 
905 #ifdef GLPAPP_SHOW_BOUNDARY_ASSEMBLY
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  //
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 
976  Epetra_SerialDenseMatrix Bt(2,2);
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 */
1110  const Epetra_IntSerialDenseVector & ipindx,
1111  const Epetra_SerialDenseMatrix & ipcoords,
1112  const Epetra_IntSerialDenseVector & pindx,
1113  const Epetra_SerialDenseMatrix & pcoords,
1114  const Epetra_IntSerialDenseMatrix & t,
1115  const Epetra_IntSerialDenseMatrix & e,
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.
1141  int format = Epetra_FECrsMatrix::COLUMN_MAJOR;
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.
1425  Epetra_SerialDenseMatrix b(2,2);
1426  Epetra_SerialDenseMatrix BtB(2,2);
1428  Epetra_SerialDenseMatrix M1(3,3);
1429  Epetra_SerialDenseMatrix M2(3,3);
1430  Epetra_SerialDenseMatrix M3(3,3);
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;
1506  Epetra_IntSerialDenseVector ipiv(dim);
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;
1532  Epetra_IntSerialDenseVector ipiv(dim);
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 
1554  Epetra_IntSerialDenseVector & ipindx,
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);
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,
1681  const Epetra_IntSerialDenseMatrix & t,
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
1708  Epetra_SerialDenseVector 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,
1845  const Epetra_IntSerialDenseMatrix & t,
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 
1863  int format = Epetra_FECrsMatrix::COLUMN_MAJOR;
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
1871  Epetra_SerialDenseVector 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,
1995  const Epetra_IntSerialDenseMatrix & t,
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
2020  Epetra_SerialDenseVector 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 
2127 bool 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 
2229 bool 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 
2297 bool 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 */
2358 int GLpApp::quadrature(const int dim, const int order,
2359  Epetra_SerialDenseMatrix & nodes,
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 }
Teuchos::RCP< Epetra_FEVector > Ny_
Teuchos::RCP< Epetra_FECrsMatrix > getR()
int NumGlobalElements() const
int GlobalLength() const
Teuchos::RCP< Epetra_FEVector > getq()
Teuchos::RCP< const Epetra_Comm > commptr_
Teuchos::RCP< Epetra_FEVector > b_
Right-hand side of the PDE.
Teuchos::RCP< const Epetra_IntSerialDenseVector > getpindx()
double beta_
Regularization parameter.
int Size(int Length_in)
void f()
void PrintVec(const Teuchos::RCP< const Epetra_Vector > &x)
Outputs the solution vector to files.
int Resize(int Length_in)
Epetra_SerialDenseMatrix NNw
int NumMyEntries(int Row) const
int Shape(int NumRows, int NumCols)
Teuchos::RCP< const Epetra_IntSerialDenseMatrix > gett()
void computeNpy(const Teuchos::RCP< const Epetra_MultiVector > &y)
Calls the function that computes the Jacobian of the nonlinear term.
Epetra_SerialDenseMatrix Nx2
int ExtractGlobalRowCopy(int_type Row, int Length, int &NumEntries, double *values, int_type *Indices) const
float NRM2(const int N, const float *X, const int INCX=1) const
Epetra_SerialDenseMatrix * NdNdx2Nw
Teuchos::RCP< Epetra_CrsMatrix > Augmat_
Augmented system matrix: [ I Jac* ] [Jac 0 ].
int MyGlobalElements(int *MyGlobalElementList) const
virtual void Print(std::ostream &os) const
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 NumGlobalRows() const
int Size(int Length_in)
Epetra_SerialDenseMatrix S3
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void computeNy(const Teuchos::RCP< const Epetra_MultiVector > &y)
Calls the function that computes the nonlinear term.
int MyLength() const
void gfctn(const Epetra_SerialDenseVector &, Epetra_SerialDenseVector &)
bool Vector2MATLAB(const Epetra_Vector &, std::ostream &)
Epetra_SerialDenseVector Weights
bool FEVector2MATLAB(const Epetra_FEVector &, std::ostream &)
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 > &)
Teuchos::RCP< Epetra_FECrsMatrix > A_
Volume stiffness matrix.
Epetra_SerialDenseMatrix S1
Teuchos::RCP< const Epetra_IntSerialDenseMatrix > gete()
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)
double determinant(const Epetra_SerialDenseMatrix &)
Provides the interface to generic abstract vector libraries.
void Print(std::ostream &os) const
bool IndicesAreLocal() const
virtual void Barrier() const =0
Teuchos::RCP< Epetra_FECrsMatrix > getNpy()
ostream & operator<<(ostream &os, const pair< Packet, Packet > &arg)
Transform to form the explicit transpose of a Epetra_RowMatrix.
Teuchos::RCP< Epetra_IntSerialDenseMatrix > t_
Elements (this includes all overlapping nodes).
void gpfctn(const Epetra_SerialDenseVector &v, Epetra_SerialDenseVector &gv)
Teuchos::RCP< Epetra_SerialDenseMatrix > pcoords_
Coordinates of all nodes in this subdomain.
virtual int MyPID() const =0
int FillComplete(bool OptimizeDataStorage=true)
int Resize(int Length_in)
int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
void GETRF(const int M, const int N, float *A, const int LDA, int *IPIV, int *INFO) const
const int * A() const
int IndexBase() const
Epetra_SerialDenseMatrix S2
virtual void Print(std::ostream &os) const
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...
double * A() const
Epetra_SerialDenseMatrix * NdNdx1Nw
int NumMyElements() const
const int N
void Print(std::ostream &os) const
Teuchos::RCP< Epetra_FECrsMatrix > Npy_
Jacobian of the nonlinear term.
Teuchos::RCP< Epetra_FECrsMatrix > R_
Edge mass matrix.
void computeAugmat()
Assembles the augmented system (KKT-type) matrix.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::RCP< const Epetra_SerialDenseMatrix > getpcoords()
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)
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 ExtractGlobalRowCopy(int GlobalRow, int Length, int &NumEntries, double *Values, int *Indices) const
T_To & dyn_cast(T_From &from)
Teuchos::RCP< Epetra_FECrsMatrix > B_
Control/state mass matrix.
void g2pfctn(const Epetra_SerialDenseVector &, Epetra_SerialDenseVector &)
virtual const Epetra_BlockMap & Map() const =0
int compproduct(Epetra_SerialDenseVector &, double *, double *)
int NumMyRows() const
const double tol
Teuchos::RCP< Epetra_SerialDenseMatrix > ipcoords_
Coordinates of nodes that are unique to this subdomain.
int GlobalAssemble(bool callFillComplete=true, Epetra_CombineMode combineMode=Add, bool save_off_and_reuse_map_exporter=false)
int lassembly(const Epetra_SerialDenseMatrix &, const Epetra_SerialDenseVector &, const Epetra_SerialDenseMatrix &, const Epetra_SerialDenseVector &, const Epetra_SerialDenseVector &, const Usr_Par &, Epetra_SerialDenseMatrix &, Epetra_SerialDenseVector &)
Teuchos::RCP< const Epetra_IntSerialDenseVector > getipindx()
Teuchos::RCP< Epetra_FECrsMatrix > getH()
void g()
int LID(int GID) const
int NumGlobalNonzeros() const
Teuchos::RCP< const Epetra_SerialDenseMatrix > getipcoords()
void computeAll(const GenSQP::Vector &x)
Calls functions to compute nonlinear quantities and the augmented system matrix.
Teuchos::RCP< Epetra_FECrsMatrix > getB()
void GETRI(const int N, float *A, const int LDA, int *IPIV, float *WORK, const int *LWORK, int *INFO) const
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 > &)
Teuchos::RCP< Epetra_FEVector > q_
The desired state.
Teuchos::RCP< Epetra_FECrsMatrix > getA()
Epetra_SerialDenseMatrix Nodes
int Scale(double ScalarA)
Teuchos::RCP< Epetra_IntSerialDenseVector > pindx_
Global nodes (interior + shared, overlapping) in this subdomain.
virtual void Print(std::ostream &os) const
int Multiply(char TransA, char TransB, double ScalarAB, const Epetra_SerialDenseMatrix &A, const Epetra_SerialDenseMatrix &B, double ScalarThis)
virtual int NumProc() const =0
int GRID(int LRID_in) const
Teuchos::RCP< Epetra_CrsMatrix > getAugmat()
bool CrsMatrix2MATLAB(const Epetra_CrsMatrix &, std::ostream &)
Epetra_SerialDenseMatrix Nx1
const Epetra_Map & DomainMap() const
Epetra_SerialDenseVector Nw
Epetra_SerialDenseMatrix * NNNw
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.
Teuchos::RCP< Epetra_FEVector > getNy()
int quadrature(const int, const int, Epetra_SerialDenseMatrix &, Epetra_SerialDenseVector &)
float DOT(const int N, const float *X, const float *Y, const int INCX=1, const int INCY=1) const
double Dot(const Epetra_SerialDenseVector &x) const
int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
std::ostream & operator<<(std::ostream &, const Usr_Par &)
The GenSQP::Vector / (y,u) Epetra_MultiVector adapter class.
double * Values() const
int Shape(int NumRows, int NumCols)
Teuchos::RCP< const Epetra_Comm > getCommPtr()
const Epetra_Comm & Comm() const
Epetra_SerialDenseMatrix N
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
Teuchos::RCP< Epetra_FEVector > getb()
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 GCID(int LCID_in) const
Teuchos::RCP< Epetra_IntSerialDenseMatrix > e_
Edges.
int inverse(const Epetra_SerialDenseMatrix &, Epetra_SerialDenseMatrix &)
Teuchos::RCP< Epetra_IntSerialDenseVector > ipindx_
Global nodes (interior, nonoverlapping) in this subdomain.
int IndexBase() const
int NumMyNonzeros() const
Teuchos::RCP< Epetra_FECrsMatrix > H_
Volume mass matrix.