GLpApp_GLpYUEpetraDataPool.cpp

Go to the documentation of this file.
00001 
00002 #include <stdlib.h>
00003 #include <algorithm>
00004 #include <iostream>
00005 
00006 #include "GLpApp_GLpYUEpetraDataPool.hpp"
00007 #include "rect2DMeshGenerator.hpp"
00008 
00009 #include "Amesos_Klu.h"
00010 #include "EpetraExt_MatrixMatrix.h"
00011 #include "EpetraExt_Reindex_LinearProblem.h"
00012 #include "EpetraExt_Transpose_RowMatrix.h"
00013 #include "Epetra_BLAS.h"
00014 #include "Epetra_CrsMatrix.h"
00015 #include "Epetra_Export.h"
00016 #include "Epetra_FECrsMatrix.h"
00017 #include "Epetra_FEVector.h"
00018 #include "Epetra_Import.h"
00019 #include "Epetra_IntSerialDenseVector.h"
00020 #include "Epetra_LAPACK.h"
00021 #include "Epetra_Map.h"
00022 #include "Epetra_MultiVector.h"
00023 #include "Epetra_SerialDenseMatrix.h"
00024 #include "Epetra_Vector.h"
00025 #include "Teuchos_ParameterList.hpp"
00026 #include "Teuchos_RefCountPtr.hpp"
00027 #include "Teuchos_TestForException.hpp"
00028 #include "Teuchos_VerboseObject.hpp"
00029 
00030 #ifdef HAVE_MPI
00031 #  include "Epetra_MpiComm.h"
00032 #  include "mpi.h"
00033 #else
00034 #  include "Epetra_SerialComm.h"
00035 #endif
00036 
00037 // Define to see all debug output for mesh generation
00038 //#define GLPYUEPETRA_DATAPOOL_DUMP_ALL_MESH
00039 
00040 namespace GLpApp {
00041 
00042 //
00043 // Declarations
00044 //
00045 
00046 const double GLp_pi = 3.14159265358979323846;
00047 
00048 ostream& operator <<(ostream &, const Usr_Par &);
00049 
00050 bool CrsMatrix2MATLAB(const Epetra_CrsMatrix &, ostream &);
00051 
00052 bool Vector2MATLAB( const Epetra_Vector &, ostream &);
00053 
00054 bool FEVector2MATLAB( const Epetra_FEVector &, ostream &);
00055 
00056 int quadrature(const int, const int, Epetra_SerialDenseMatrix &,
00057                Epetra_SerialDenseVector &);
00058 
00059 int meshreader(
00060   const Epetra_Comm &,
00061   Epetra_IntSerialDenseVector &,
00062   Epetra_SerialDenseMatrix &,
00063   Epetra_IntSerialDenseVector &,
00064   Epetra_SerialDenseMatrix &,
00065   Epetra_IntSerialDenseMatrix &,
00066   Epetra_IntSerialDenseMatrix &,
00067   const char geomFileBase[],
00068   const bool trace = true,
00069   const bool dumpAll = false
00070   );
00071 
00072 int lassembly(const Epetra_SerialDenseMatrix &,
00073               const Epetra_SerialDenseVector &,
00074               const Epetra_SerialDenseMatrix &,
00075               const Epetra_SerialDenseVector &,
00076               const Epetra_SerialDenseVector &,
00077               const Usr_Par &,
00078               Epetra_SerialDenseMatrix &,
00079               Epetra_SerialDenseVector &);
00080 
00081 int assemblyFECrs(const Epetra_Comm &,
00082                   const Epetra_IntSerialDenseVector &,
00083                   const Epetra_SerialDenseMatrix &,
00084                   const Epetra_IntSerialDenseVector &,
00085                   const Epetra_SerialDenseMatrix &,
00086                   const Epetra_IntSerialDenseMatrix &,
00087                   const Epetra_IntSerialDenseMatrix &,
00088                   Teuchos::RefCountPtr<Epetra_FECrsMatrix> &,
00089                   Teuchos::RefCountPtr<Epetra_FEVector> &);
00090 
00091 int assemblyFECrs(const Epetra_Comm &,
00092                   const Epetra_IntSerialDenseVector &,
00093                   const Epetra_SerialDenseMatrix &,
00094                   const Epetra_IntSerialDenseVector &,
00095                   const Epetra_SerialDenseMatrix &,
00096                   const Epetra_IntSerialDenseMatrix &,
00097                   const Epetra_IntSerialDenseMatrix &,
00098                   Teuchos::RefCountPtr<Epetra_FECrsMatrix> &,
00099                   Teuchos::RefCountPtr<Epetra_FEVector> &,
00100                   bool);
00101 
00102 int assemble(const Epetra_Comm &,
00103              const Epetra_IntSerialDenseVector &,
00104              const Epetra_SerialDenseMatrix &,
00105              const Epetra_IntSerialDenseVector &,
00106              const Epetra_SerialDenseMatrix &,
00107              const Epetra_IntSerialDenseMatrix &,
00108              const Epetra_IntSerialDenseMatrix &,
00109              Teuchos::RefCountPtr<Epetra_FECrsMatrix> &,
00110              Teuchos::RefCountPtr<Epetra_FECrsMatrix> &,
00111              Teuchos::RefCountPtr<Epetra_FEVector> &);
00112 
00113 int assemble_bdry(
00114   const Epetra_Comm                                &Comm
00115   ,const Epetra_IntSerialDenseVector               &ipindx
00116   ,const Epetra_SerialDenseMatrix                  &ipcoords
00117   ,const Epetra_IntSerialDenseVector               &pindx
00118   ,const Epetra_SerialDenseMatrix                  &pcoords
00119   ,const Epetra_IntSerialDenseMatrix               &t
00120   ,const Epetra_IntSerialDenseMatrix               &e
00121   ,Teuchos::RefCountPtr<Epetra_FECrsMatrix>        *B
00122   ,Teuchos::RefCountPtr<Epetra_FECrsMatrix>        *R
00123   );
00124 
00125 int nonlinvec(const Epetra_Comm &,
00126               const Epetra_IntSerialDenseVector &,
00127               const Epetra_SerialDenseMatrix &,
00128               const Epetra_IntSerialDenseVector &,
00129               const Epetra_SerialDenseMatrix &,
00130               const Epetra_IntSerialDenseMatrix &,
00131               const Teuchos::RefCountPtr<const Epetra_MultiVector> &,
00132               Teuchos::RefCountPtr<Epetra_FEVector> &);
00133 
00134 int nonlinjac(const Epetra_Comm &,
00135               const Epetra_IntSerialDenseVector &,
00136               const Epetra_SerialDenseMatrix &,
00137               const Epetra_IntSerialDenseVector &,
00138               const Epetra_SerialDenseMatrix &,
00139               const Epetra_IntSerialDenseMatrix &,
00140               const Teuchos::RefCountPtr<const Epetra_MultiVector> &,
00141               Teuchos::RefCountPtr<Epetra_FECrsMatrix> &);
00142 
00143 int nonlinhessvec(const Epetra_Comm &,
00144                   const Epetra_IntSerialDenseVector &,
00145                   const Epetra_SerialDenseMatrix &,
00146                   const Epetra_IntSerialDenseVector &,
00147                   const Epetra_SerialDenseMatrix &,
00148                   const Epetra_IntSerialDenseMatrix &,
00149                   const Teuchos::RefCountPtr<const Epetra_MultiVector> &,
00150                   const Teuchos::RefCountPtr<const Epetra_MultiVector> &,
00151                   const Teuchos::RefCountPtr<const Epetra_MultiVector> &,
00152                   Teuchos::RefCountPtr<Epetra_FEVector> &);
00153 
00154 int compproduct(Epetra_SerialDenseVector &, double *, double *);
00155 
00156 int compproduct(Epetra_SerialDenseVector &, double *, double *, double *);
00157 
00158 double determinant(const Epetra_SerialDenseMatrix &);
00159 
00160 int inverse(const Epetra_SerialDenseMatrix &, Epetra_SerialDenseMatrix &);
00161 
00162 int quadrature(
00163   const int, const int, Epetra_SerialDenseMatrix &,
00164   Epetra_SerialDenseVector &);
00165 
00166 void gpfctn(const Epetra_SerialDenseVector & v, Epetra_SerialDenseVector & gv);
00167 
00168 void g2pfctn(const Epetra_SerialDenseVector & , Epetra_SerialDenseVector & );
00169 
00170 void gfctn(const Epetra_SerialDenseVector & , Epetra_SerialDenseVector & );
00171 
00172 //
00173 // GLpYUEpetraDataPool
00174 //
00175 
00176 GLpYUEpetraDataPool::GLpYUEpetraDataPool(
00177   Teuchos::RefCountPtr<const Epetra_Comm>    const& commptr
00178   ,const double                              beta
00179   ,const double                              len_x
00180   ,const double                              len_y
00181   ,const int                                 local_nx
00182   ,const int                                 local_ny
00183   ,const char                                myfile[]
00184   ,const bool                                trace
00185   )
00186   :commptr_(commptr)
00187   ,beta_(beta)
00188 {
00189   ipcoords_ = Teuchos::rcp( new Epetra_SerialDenseMatrix() );
00190   ipindx_ = Teuchos::rcp( new Epetra_IntSerialDenseVector() );
00191   pcoords_ = Teuchos::rcp( new Epetra_SerialDenseMatrix() );
00192   pindx_ = Teuchos::rcp( new Epetra_IntSerialDenseVector() );
00193   t_ = Teuchos::rcp( new Epetra_IntSerialDenseMatrix() );
00194   e_ = Teuchos::rcp( new Epetra_IntSerialDenseMatrix() );
00195 
00196   if( myfile && myfile[0] != '\0' ) {
00197     meshreader(*commptr_,*ipindx_,*ipcoords_,*pindx_,*pcoords_,*t_,*e_,myfile,trace);
00198   }
00199   else {
00200     rect2DMeshGenerator(
00201       commptr_->NumProc(),commptr_->MyPID()
00202       ,len_x,len_y,local_nx,local_ny,2 // 2==Neuman boundary conditions!
00203       ,&*ipindx_,&*ipcoords_,&*pindx_,&*pcoords_,&*t_,&*e_
00204 #ifdef GLPYUEPETRA_DATAPOOL_DUMP_ALL_MESH
00205       ,&*Teuchos::VerboseObjectBase::getDefaultOStream(),true
00206 #else
00207       ,NULL,false
00208 #endif
00209       );
00210   }
00211 
00212   // Assemble volume and boundary mass and stiffness matrices, and the right-hand side of the PDE.
00213   assemble( *commptr, *ipindx_, *ipcoords_, *pindx_, *pcoords_, *t_, *e_, A_, H_, b_ );
00214   assemble_bdry( *commptr, *ipindx_, *ipcoords_, *pindx_, *pcoords_, *t_, *e_, &B_, &R_ );
00215 
00216   // Set desired state q.
00217   Epetra_Map standardmap(A_->DomainMap());
00218   q_ = Teuchos::rcp(new Epetra_FEVector(standardmap));
00219   int * qintvalues = new int[standardmap.NumMyElements()];
00220   double * qdvalues = new double[standardmap.NumMyElements()];
00221   standardmap.MyGlobalElements(qintvalues);
00222   for (int i = 0; i < standardmap.NumMyElements(); i++)
00223       qdvalues[i]=cos( GLp_pi* ((*ipcoords_)(i,0)) ) * cos( GLp_pi* ((*ipcoords_)(i,1)) );
00224   q_->ReplaceGlobalValues(standardmap.NumMyElements(), qintvalues, qdvalues);
00225   q_->GlobalAssemble();
00226 }
00227 
00228 void GLpYUEpetraDataPool::computeAll( const GenSQP::Vector &x )
00229 {
00230 
00231   // Dynamic cast back to Epetra vectors here.
00232   Teuchos::RefCountPtr<const Epetra_MultiVector> ey =
00233         (Teuchos::dyn_cast<GenSQP::YUEpetraVector>(const_cast<GenSQP::Vector&>(x))).getYVector();
00234 
00235   computeNy(ey);
00236 
00237   computeNpy(ey);
00238 
00239   computeAugmat();
00240   
00241 }
00242 
00243 int GLpYUEpetraDataPool::solveAugsys(
00244   const Teuchos::RefCountPtr<const Epetra_MultiVector> & rhsy,
00245   const Teuchos::RefCountPtr<const Epetra_MultiVector> & rhsu,
00246   const Teuchos::RefCountPtr<const Epetra_MultiVector> & rhsp,
00247   const Teuchos::RefCountPtr<Epetra_MultiVector> & y,
00248   const Teuchos::RefCountPtr<Epetra_MultiVector> & u,
00249   const Teuchos::RefCountPtr<Epetra_MultiVector> & p,
00250   double tol )
00251 {
00252 /*
00253   int systemChoice = 1;   // 1 for full KKT system solve, 2 for Schur complement solve
00254   int solverChoice = 12;  // These options are for the full KKT system solve.
00255                           // 11 for AztecOO with built-in Schwarz DD preconditioning and ILU on subdomains
00256                           // 12 for AztecOO with IFPACK Schwarz DD preconditioning and Umfpack on subdomains
00257                           // 13 for a direct sparse solver (Umfpack, KLU)
00258   
00259   if (systemChoice == 1) {
00260     // We're using the full KKT system formulation to solve the augmented system.
00261    
00262     Epetra_Map standardmap(A_->DomainMap());
00263     int numstates = standardmap.NumGlobalElements();
00264     Epetra_Map bdryctrlmap(B_->DomainMap());
00265     int numcontrols = bdryctrlmap.NumGlobalElements();
00266     Epetra_Vector rhs( (Epetra_BlockMap&)Augmat_->RangeMap() );
00267     Epetra_Vector soln( (Epetra_BlockMap&)Augmat_->RangeMap() );
00268     soln.Random();  
00269 
00270     std::vector<double> values(rhsy->MyLength() + rhsu->MyLength() + rhsp->MyLength());
00271     std::vector<int> indices(rhsy->MyLength() + rhsu->MyLength() + rhsp->MyLength());
00272     ((Epetra_BlockMap&)Augmat_->RangeMap()).MyGlobalElements(&indices[0]);
00273 
00274     for (int i=0; i<rhsy->MyLength(); i++) {
00275       values[i] = (*((*rhsy)(0)))[i];
00276     }
00277     for (int i=0; i<rhsu->MyLength(); i++) {
00278       values[i+rhsy->MyLength()] = (*((*rhsu)(0)))[i];
00279     }
00280     for (int i=0; i<rhsp->MyLength(); i++) {
00281       values[i+rhsy->MyLength()+rhsu->MyLength()] = (*((*rhsp)(0)))[i];
00282     }
00283 
00284     rhs.ReplaceGlobalValues(rhsy->MyLength() + rhsu->MyLength() + rhsp->MyLength(), &values[0], &indices[0]);
00285 
00286     if (solverChoice == 11) {
00287       int Overlap = 3;
00288       int ival = 4;
00289 
00290       AztecOO::AztecOO kktsolver(&(*Augmat_), &soln, &rhs);
00291       kktsolver.SetAztecOption( AZ_solver, AZ_gmres );
00292       kktsolver.SetAztecOption( AZ_precond, AZ_dom_decomp );
00293       kktsolver.SetAztecOption(AZ_subdomain_solve, AZ_ilu);
00294       //kktsolver.SetAztecOption( AZ_kspace, 2*numstates+numcontrols );
00295       kktsolver.SetAztecOption( AZ_kspace, 9000 );
00296       kktsolver.SetAztecOption(AZ_overlap,Overlap);
00297       kktsolver.SetAztecOption(AZ_graph_fill,ival);
00298       //kktsolver.SetAztecOption(AZ_poly_ord, ival);
00299       //kktsolver.SetAztecParam(AZ_drop, 1e-9);
00300       kktsolver.SetAztecParam(AZ_athresh, 1e-5);
00301       //kktsolver.SetAztecParam(AZ_rthresh, 0.0);
00302       kktsolver.SetAztecOption( AZ_reorder, 0 );
00303       //kktsolver.SetAztecParam44( AZ_ilut_fill, 1.5 );
00304       kktsolver.SetAztecOption( AZ_output, AZ_last );
00305       //kktsolver.Iterate(2*numstates+numcontrols,1e-12);
00306       kktsolver.Iterate(9000,1e-11);
00307       //cout << soln;
00308     }
00309     else if (solverChoice == 12) {
00310       // =============================================================== //
00311       // 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 //
00312       // =============================================================== //
00313 
00314       Teuchos::ParameterList List;
00315 
00316       // allocates an IFPACK factory. No data is associated
00317       // to this object (only method Create()).
00318       Ifpack Factory;
00319 
00320       // create the preconditioner. For valid PrecType values,
00321       // please check the documentation
00322       string PrecType = "Amesos";
00323       int OverlapLevel = 2; // must be >= 0. If Comm.NumProc() == 1,
00324                             // it is ignored.
00325   
00326       Ifpack_Preconditioner* Prec = Factory.Create(PrecType, &(*Augmat_), OverlapLevel);
00327       assert(Prec != 0);
00328 
00329       // specify the Amesos solver to be used.
00330       // If the selected solver is not available,
00331       // IFPACK will try to use Amesos' KLU (which is usually always
00332       // compiled). Amesos' serial solvers are:
00333       // "Amesos_Klu", "Amesos_Umfpack", "Amesos_Superlu"
00334       List.set("amesos: solver type", "Amesos_Umfpack");
00335 
00336       // sets the parameters
00337       IFPACK_CHK_ERR(Prec->SetParameters(List));
00338 
00339       // initialize the preconditioner. At this point the matrix must
00340       // have been FillComplete()'d, but actual values are ignored.
00341       // At this call, Amesos will perform the symbolic factorization.
00342       IFPACK_CHK_ERR(Prec->Initialize());
00343 
00344       // Builds the preconditioners, by looking for the values of
00345       // the matrix. At this call, Amesos will perform the
00346       // numeric factorization.
00347       IFPACK_CHK_ERR(Prec->Compute());
00348 
00349       // =================================================== //
00350       // E N D   O F   I F P A C K   C O N S T R U C T I O N //
00351       // =================================================== //
00352 
00353       // need an Epetra_LinearProblem to define AztecOO solver
00354       Epetra_LinearProblem Problem;
00355       Problem.SetOperator(&(*Augmat_));
00356       Problem.SetLHS(&soln);
00357       Problem.SetRHS(&rhs);
00358 
00359       // now we can allocate the AztecOO solver
00360       AztecOO kktsolver(Problem);
00361 
00362       // specify solver
00363       kktsolver.SetAztecOption(AZ_solver,AZ_gmres);
00364       kktsolver.SetAztecOption(AZ_kspace, 300 );
00365       kktsolver.SetAztecOption(AZ_output,AZ_last);
00366 
00367       // HERE WE SET THE IFPACK PRECONDITIONER
00368       kktsolver.SetPrecOperator(Prec);
00369 
00370       // .. and here we solve
00371       kktsolver.Iterate(300,1e-12);
00372 
00373       // delete the preconditioner
00374       delete Prec;
00375     }
00376     else if (solverChoice == 13) {
00377       Epetra_LinearProblem Problem;
00378       Problem.SetOperator(&(*Augmat_));
00379       Problem.SetLHS(&soln);
00380       Problem.SetRHS(&rhs);
00381       
00382       EpetraExt::LinearProblem_Reindex reindex(NULL);
00383       Epetra_LinearProblem newProblem = reindex(Problem);
00384       
00385       Amesos_Klu kktsolver(newProblem);
00386    
00387       AMESOS_CHK_ERR(kktsolver.SymbolicFactorization());
00388       AMESOS_CHK_ERR(kktsolver.NumericFactorization());
00389       AMESOS_CHK_ERR(kktsolver.Solve());
00390       kktsolver.PrintTiming();
00391     }
00392     
00393     
00394     for (int i=0; i<rhsy->MyLength(); i++) {
00395       (*((*y)(0)))[i] = soln[i];
00396     }
00397     for (int i=0; i<rhsu->MyLength(); i++) {
00398       (*((*u)(0)))[i] = soln[i+rhsy->MyLength()];
00399     }
00400     for (int i=0; i<rhsp->MyLength(); i++) {
00401       (*((*p)(0)))[i] = soln[i+rhsy->MyLength()+rhsu->MyLength()];
00402     }
00403     
00404   }
00405   else if (systemChoice == 2) {
00406     // We're using the Schur complement formulation to solve the augmented system.
00407   
00408     // Form linear operator.
00409     GLpApp::SchurOp schurop(A_, B_, Npy_);
00410   
00411     // Form Schur complement right-hand side.
00412     Epetra_MultiVector ny( (Epetra_BlockMap&)Npy_->RangeMap(), 1);
00413     Epetra_MultiVector ay( (Epetra_BlockMap&)A_->RangeMap(), 1);
00414     Epetra_MultiVector schurrhs( (Epetra_BlockMap&)A_->RangeMap(), 1);
00415     Epetra_MultiVector bu( (Epetra_BlockMap&)B_->RangeMap(), 1);
00416     A_->Multiply(false, *rhsy, ay);
00417     Npy_->Multiply(false, *rhsy, ny);
00418     B_->Multiply(false, *rhsu, bu);
00419     schurrhs.Update(1.0, ny, 1.0, ay, 0.0);
00420     schurrhs.Update(-1.0, *rhsp, 1.0, bu, 1.0);
00421   
00422     p->PutScalar(0.0);
00423     Epetra_LinearProblem linprob(&schurop, &(*p), &schurrhs);
00424     AztecOO::AztecOO Solver(linprob);
00425     Solver.SetAztecOption( AZ_solver, AZ_cg );
00426     Solver.SetAztecOption( AZ_precond, AZ_none );
00427     Solver.SetAztecOption( AZ_output, AZ_none );
00428     Solver.Iterate(8000,tol);
00429   
00430     Epetra_MultiVector bp( (Epetra_BlockMap&)B_->DomainMap(), 1);
00431     B_->Multiply(true, *p, bp);
00432     u->Update(1.0, *rhsu, -1.0, bp, 0.0);
00433 
00434     Epetra_MultiVector ap( (Epetra_BlockMap&)A_->DomainMap(), 1);
00435     Epetra_MultiVector np( (Epetra_BlockMap&)A_->DomainMap(), 1);
00436     A_->Multiply(true, *p, ap);
00437     Npy_->Multiply(true, *p, np);
00438     y->Update(1.0, *rhsy,0.0);
00439     y->Update(-1.0, ap, -1.0, np, 1.0);
00440   }
00441 */
00442   TEST_FOR_EXCEPT(true);
00443   return 0;
00444 }
00445 
00446 Teuchos::RefCountPtr<const Epetra_Comm> GLpYUEpetraDataPool::getCommPtr()   { return commptr_; }
00447 
00448 Teuchos::RefCountPtr<Epetra_FECrsMatrix> GLpYUEpetraDataPool::getA()  { return A_; }
00449 
00450 Teuchos::RefCountPtr<Epetra_FECrsMatrix> GLpYUEpetraDataPool::getB()  { return B_; }
00451 
00452 Teuchos::RefCountPtr<Epetra_FECrsMatrix> GLpYUEpetraDataPool::getH()  { return H_; }
00453 
00454 Teuchos::RefCountPtr<Epetra_FECrsMatrix> GLpYUEpetraDataPool::getR()  { return R_; }
00455 
00456 Teuchos::RefCountPtr<Epetra_CrsMatrix> GLpYUEpetraDataPool::getAugmat()  { return Augmat_; }
00457 
00458 Teuchos::RefCountPtr<Epetra_FECrsMatrix> GLpYUEpetraDataPool::getNpy()  { return Npy_; }
00459 
00460 Teuchos::RefCountPtr<Epetra_FEVector> GLpYUEpetraDataPool::getb()  { return b_; }
00461 
00462 Teuchos::RefCountPtr<Epetra_FEVector> GLpYUEpetraDataPool::getq()  { return q_; }
00463 
00464 Teuchos::RefCountPtr<Epetra_FEVector> GLpYUEpetraDataPool::getNy()  { return Ny_; }
00465 
00466 double GLpYUEpetraDataPool::getbeta()  { return beta_; }
00467 
00468 Teuchos::RefCountPtr<const Epetra_SerialDenseMatrix> GLpYUEpetraDataPool::getipcoords()  { return ipcoords_; }
00469 
00470 Teuchos::RefCountPtr<const Epetra_IntSerialDenseVector> GLpYUEpetraDataPool::getipindx()  { return ipindx_; }
00471 
00472 Teuchos::RefCountPtr<const Epetra_SerialDenseMatrix> GLpYUEpetraDataPool::getpcoords()  { return pcoords_; }
00473 
00474 Teuchos::RefCountPtr<const Epetra_IntSerialDenseVector> GLpYUEpetraDataPool::getpindx()  { return pindx_; }
00475 
00476 Teuchos::RefCountPtr<const Epetra_IntSerialDenseMatrix> GLpYUEpetraDataPool::gett()  { return t_; }
00477 
00478 Teuchos::RefCountPtr<const Epetra_IntSerialDenseMatrix> GLpYUEpetraDataPool::gete()  { return e_; }
00479 
00480 
00481 void GLpYUEpetraDataPool::computeNy( const Teuchos::RefCountPtr<const Epetra_MultiVector> & y )
00482 {
00483   Epetra_Map overlapmap(-1, pindx_->M(), (int*)(pindx_)->A(), 1, *commptr_);
00484   Epetra_Map standardmap(A_->DomainMap());
00485   Teuchos::RefCountPtr<Epetra_MultiVector> yoverlap = Teuchos::rcp(new Epetra_MultiVector(overlapmap, 1));
00486   Epetra_Import Importer(overlapmap, standardmap);
00487   yoverlap->Import(*y, Importer, Insert);
00488   nonlinvec(*commptr_, *ipindx_, *ipcoords_, *pindx_, *pcoords_, *t_, yoverlap, Ny_);
00489 }
00490 
00491 
00492 void GLpYUEpetraDataPool::computeNpy( const Teuchos::RefCountPtr<const Epetra_MultiVector> & y )
00493 {
00494   Epetra_Map overlapmap(-1, pindx_->M(), (int*)(pindx_)->A(), 1, *commptr_);
00495   Epetra_Map standardmap(A_->DomainMap());
00496   Teuchos::RefCountPtr<Epetra_MultiVector> yoverlap = Teuchos::rcp(new Epetra_MultiVector(overlapmap, 1));
00497   Epetra_Import Importer(overlapmap, standardmap);
00498   yoverlap->Import(*y, Importer, Insert);
00499   nonlinjac(*commptr_, *ipindx_, *ipcoords_, *pindx_, *pcoords_, *t_, yoverlap, Npy_);
00500 }
00501 
00502 
00503 void GLpYUEpetraDataPool::computeAugmat()
00504 {
00505   Epetra_Map standardmap(A_->DomainMap());
00506   Epetra_Map bdryctrlmap(B_->DomainMap());
00507 
00508   int indexBase = 1;
00509 
00510   int numstates = standardmap.NumGlobalElements();
00511   //int numcontrols = bdryctrlmap.NumGlobalElements();
00512   int nummystates = standardmap.NumMyElements();
00513   int nummycontrols = bdryctrlmap.NumMyElements();
00514 
00515   Epetra_IntSerialDenseVector KKTmapindx(2*nummystates+nummycontrols);
00516   
00517   // Build KKT map.
00518   Epetra_IntSerialDenseVector states(nummystates);
00519   Epetra_IntSerialDenseVector controls(nummycontrols);
00520   standardmap.MyGlobalElements(states.Values());
00521   bdryctrlmap.MyGlobalElements(controls.Values());
00522   for (int i=0; i<nummystates; i++) {
00523     KKTmapindx(i) = states(i);
00524     KKTmapindx(nummystates+nummycontrols+i) = 2*numstates+states(i);
00525   }
00526   for (int i=0; i<nummycontrols; i++) {
00527     KKTmapindx(nummystates+i) = numstates+controls(i);
00528   }
00529   Epetra_Map KKTmap(-1, 2*nummystates+nummycontrols, KKTmapindx.Values(), indexBase, *(commptr_));
00530   
00531   
00532   // Start building the KKT matrix.
00533   
00534   Augmat_ = Teuchos::rcp(new Epetra_CrsMatrix(Copy, KKTmap, 0));
00535 
00536   double one[1];
00537   one[0]=1.0;
00538   for (int i=0; i<nummystates+nummycontrols; i++) {
00539     Augmat_->InsertGlobalValues(KKTmapindx.Values()[i], 1, one, KKTmapindx.Values()+i);
00540   }
00541   
00542   int checkentries=0;
00543   int nummyentries=0;
00544   Epetra_SerialDenseVector values(nummyentries);
00545   Epetra_IntSerialDenseVector indices(nummyentries);
00546   // Insert A and Npy into Augmat.
00547   for (int i=0; i<nummystates; i++) {
00548     nummyentries = A_->NumMyEntries(i);
00549     values.Resize(nummyentries);
00550     indices.Resize(nummyentries);
00551     A_->ExtractGlobalRowCopy(KKTmapindx.Values()[i], nummyentries, checkentries, values.Values(),
00552                              indices.Values());
00553     if (nummyentries > 0)
00554       Augmat_->InsertGlobalValues(KKTmapindx.Values()[i]+2*numstates, nummyentries, values.Values(), 
00555                                   indices.Values());
00556     nummyentries = Npy_->NumMyEntries(i);
00557     values.Resize(nummyentries);
00558     indices.Resize(nummyentries);
00559     Npy_->ExtractGlobalRowCopy(KKTmapindx.Values()[i], nummyentries, checkentries, values.Values(),
00560                              indices.Values());
00561     if (nummyentries > 0)
00562       Augmat_->InsertGlobalValues(KKTmapindx.Values()[i]+2*numstates, nummyentries, values.Values(), 
00563                                   indices.Values());
00564   }
00565   // Insert B into Augmat.
00566   for (int i=0; i<nummystates; i++) {
00567     nummyentries = B_->NumMyEntries(i);
00568     values.Resize(nummyentries);
00569     indices.Resize(nummyentries);
00570     B_->ExtractGlobalRowCopy(KKTmapindx.Values()[i], nummyentries, checkentries, values.Values(),
00571                              indices.Values());
00572     for (int j=0; j<nummyentries; j++)
00573       indices[j] = indices[j]+numstates;
00574     if (nummyentries > 0)
00575       Augmat_->InsertGlobalValues(KKTmapindx.Values()[i]+2*numstates, nummyentries, values.Values(), 
00576                                   indices.Values());
00577   }
00578   
00579   bool MakeDataContiguous = false;
00580   EpetraExt::RowMatrix_Transpose transposer( MakeDataContiguous );
00581   Epetra_CrsMatrix & transA = dynamic_cast<Epetra_CrsMatrix&>(transposer(*A_));
00582   Epetra_CrsMatrix & transB = dynamic_cast<Epetra_CrsMatrix&>(transposer(*B_));
00583   Epetra_CrsMatrix & transNpy = dynamic_cast<Epetra_CrsMatrix&>(transposer(*Npy_));
00584   // Insert transpose of A and Npy into Augmat.
00585   for (int i=0; i<nummystates; i++) {
00586     nummyentries = transA.NumMyEntries(i);
00587     values.Resize(nummyentries);
00588     indices.Resize(nummyentries);
00589     transA.ExtractGlobalRowCopy(KKTmapindx.Values()[i], nummyentries, checkentries, values.Values(),
00590                                 indices.Values());
00591     for (int j=0; j<nummyentries; j++)
00592       indices[j] = indices[j]+2*numstates;
00593     if (nummyentries > 0)
00594       Augmat_->InsertGlobalValues(KKTmapindx.Values()[i], nummyentries, values.Values(), 
00595                                   indices.Values());
00596     nummyentries = transNpy.NumMyEntries(i);
00597     values.Resize(nummyentries);
00598     indices.Resize(nummyentries);
00599     transNpy.ExtractGlobalRowCopy(KKTmapindx.Values()[i], nummyentries, checkentries, values.Values(),
00600                                   indices.Values());
00601     for (int j=0; j<nummyentries; j++)
00602       indices[j] = indices[j]+2*numstates;
00603     if (nummyentries > 0)
00604       Augmat_->InsertGlobalValues(KKTmapindx.Values()[i], nummyentries, values.Values(), 
00605                                   indices.Values());
00606   }
00607   // Insert transpose of B into Augmat.
00608   for (int i=0; i<nummystates; i++) {
00609     nummyentries = transB.NumMyEntries(i);
00610     values.Resize(nummyentries);
00611     indices.Resize(nummyentries);
00612     transB.ExtractGlobalRowCopy(KKTmapindx.Values()[i], nummyentries, checkentries, values.Values(),
00613                                 indices.Values());
00614     for (int j=0; j<nummyentries; j++)
00615       indices[j] = indices[j]+2*numstates;
00616     int err = 0;
00617     if (nummyentries > 0)
00618       err = Augmat_->InsertGlobalValues(KKTmapindx.Values()[i]+numstates, nummyentries,
00619                                         values.Values(), indices.Values());
00620     // This will give a nasty message if something goes wrong with the insertion of B transpose.
00621     if (err < 0) {
00622       cout << "Insertion of entries failed:\n";
00623       cout << indices;
00624       cout << nummyentries << endl;
00625       cout << "at row: " << KKTmapindx.Values()[i]+numstates << endl << endl;
00626     }
00627   }
00628 
00629   Augmat_->FillComplete(KKTmap, KKTmap);
00630   // End building the KKT matrix.
00631 
00632 }
00633 
00634 void GLpYUEpetraDataPool::PrintVec( const Teuchos::RefCountPtr<const Epetra_Vector> & x )
00635 {
00636   Vector2MATLAB(*x, cout);
00637 }
00638 
00639 Usr_Par::Usr_Par() {
00640   Epetra_BLAS blas;
00641   Epetra_SerialDenseVector product(4);
00642 
00643   // get nodes and weights
00644   quadrature(2,3,Nodes,Weights);
00645     
00646   // Evaluate nodal basis functions and their derivatives at quadrature
00647   // pts N(i,j) = value of the j-th basis function at quadrature node i.
00648   N.Shape(Nodes.M(),3);
00649   for (int i=0; i<Nodes.M(); i++) {
00650     N(i,0) = 1.0 - Nodes(i,0) - Nodes(i,1);
00651     N(i,1) = Nodes(i,0);
00652     N(i,2) = Nodes(i,1);
00653   }
00654   // Nx1(i,j) partial derrivative of basis function j wrt x1 at node i
00655   Nx1.Shape(Nodes.M(),3);
00656   for (int i=0; i<Nodes.M(); i++) {
00657     Nx1(i,0) = -1.0;
00658     Nx1(i,1) = 1.0;
00659     Nx1(i,2) = 0.0;
00660   }
00661   // Nx2(i,j) partial derrivative of basis function j wrt x2 at node i
00662   Nx2.Shape(Nodes.M(),3);
00663   for (int i=0; i<Nodes.M(); i++) {
00664     Nx2(i,0) = -1.0;
00665     Nx2(i,1) = 0.0;
00666     Nx2(i,2) = 1.0;
00667   }
00668 
00669   // Evaluate integrals of various combinations of partial derivatives
00670   // of the local basis functions (they're constant).
00671   S1.Shape(3,3);
00672   S1(0,0)= 1.0; S1(0,1)=-1.0; S1(0,2)= 0.0;
00673   S1(1,0)=-1.0; S1(1,1)= 1.0; S1(1,2)= 0.0;
00674   S1(2,0)= 0.0; S1(2,1)= 0.0; S1(2,2)= 0.0;
00675   S2.Shape(3,3);
00676   S2(0,0)= 2.0; S2(0,1)=-1.0; S2(0,2)=-1.0;
00677   S2(1,0)=-1.0; S2(1,1)= 0.0; S2(1,2)= 1.0;
00678   S2(2,0)=-1.0; S2(2,1)= 1.0; S2(2,2)= 0.0;
00679   S3.Shape(3,3);
00680   S3(0,0)= 1.0; S3(0,1)= 0.0; S3(0,2)=-1.0;
00681   S3(1,0)= 0.0; S3(1,1)= 0.0; S3(1,2)= 0.0;
00682   S3(2,0)=-1.0; S3(2,1)= 0.0; S3(2,2)= 1.0;
00683     
00684   // Evaluate integrals of basis functions N(i).
00685   Nw.Size(3);
00686   Nw.Multiply('T', 'N', 1.0, N, Weights, 0.0);
00687 
00688   // Evaluate integrals of 9 products N(i)*N(j).
00689   NNw.Shape(3,3);
00690   for (int i=0; i<3; i++) {
00691     for (int j=0; j<3; j++) {
00692       compproduct(product, N[i], N[j]);
00693       NNw(i,j) = blas.DOT(Weights.M(), Weights.A(), product.A());
00694     }
00695   }
00696 
00697   // Evaluate integrals of 27 products N(i)*N(j)*N(k).
00698   NNNw = new Epetra_SerialDenseMatrix[3];
00699   NNNw[0].Shape(3,3); NNNw[1].Shape(3,3); NNNw[2].Shape(3,3); 
00700   for (int i=0; i<3; i++) {
00701     for (int j=0; j<3; j++) {
00702       for (int k=0; k<3; k++) {
00703         compproduct(product, N[i], N[j], N[k]);
00704         NNNw[k](i,j) = blas.DOT(Weights.M(), Weights.A(), product.A());
00705       }
00706     }
00707   }
00708 
00709   // Evaluate integrals of 27 products N(i)*(dN(j)/dx1)*N(k).
00710   NdNdx1Nw = new Epetra_SerialDenseMatrix[3];
00711   NdNdx1Nw[0].Shape(3,3); NdNdx1Nw[1].Shape(3,3); NdNdx1Nw[2].Shape(3,3); 
00712   for (int i=0; i<3; i++) {
00713     for (int j=0; j<3; j++) {
00714       for (int k=0; k<3; k++) {
00715         compproduct(product, N[i], Nx1[j], N[k]);
00716         NdNdx1Nw[k](i,j) = blas.DOT(Weights.M(), Weights.A(), product.A());
00717       }
00718     }
00719   }
00720 
00721   // Evaluate integrals of 27 products N(i)*(dN(j)/dx2)*N(k).
00722   NdNdx2Nw = new Epetra_SerialDenseMatrix[3];
00723   NdNdx2Nw[0].Shape(3,3); NdNdx2Nw[1].Shape(3,3); NdNdx2Nw[2].Shape(3,3); 
00724   for (int i=0; i<3; i++) {
00725     for (int j=0; j<3; j++) {
00726       for (int k=0; k<3; k++) {
00727         compproduct(product, N[i], Nx2[j], N[k]);
00728         NdNdx2Nw[k](i,j) = blas.DOT(Weights.M(), Weights.A(), product.A());
00729       }
00730     }
00731   }
00732 
00733 }
00734 
00735 void Usr_Par::Print(ostream& os) const {
00736   os << endl;
00737   os << "\n\nQuadrature nodes:";
00738   os << "\n-----------------";
00739   Nodes.Print(os);
00740   os << "\n\nQuadrature weights:";
00741   os << "\n-------------------\n";
00742   Weights.Print(os);
00743   os << "\n\nNodal basis functions:";
00744   os << "\n----------------------";
00745   N.Print(os);
00746   os << "\n\nIntegrals of combinations of partial derivatives:";
00747   os << "\n-------------------------------------------------";
00748   S1.Print(os);
00749   S2.Print(os);
00750   S3.Print(os);
00751   os << "\n\nIntegrals of basis functions:";
00752   os << "\n-----------------------------\n";
00753   Nw.Print(os);
00754   os << "\n\nIntegrals of products N(i)*N(j):";
00755   os << "\n--------------------------------\n";
00756   NNw.Print(os);
00757   os << "\n\nIntegrals of products N(i)*N(j)*N(k):";
00758   os << "\n-------------------------------------\n";
00759   NNNw[0].Print(os); NNNw[1].Print(os); NNNw[2].Print(os);
00760   os << "\n\nIntegrals of products N(i)*(dN(j)/dx1)*N(k):";
00761   os << "\n--------------------------------------------\n";
00762   NdNdx1Nw[0].Print(os); NdNdx1Nw[1].Print(os); NdNdx1Nw[2].Print(os);
00763   os << "\n\nIntegrals of products N(i)*(dN(j)/dx2)*N(k):";
00764   os << "\n--------------------------------------------\n";
00765   NdNdx2Nw[0].Print(os); NdNdx2Nw[1].Print(os); NdNdx2Nw[2].Print(os);
00766 }
00767 
00768 ostream& operator <<(ostream & out, const Usr_Par & usr_par)
00769 {
00770   usr_par.Print(out);
00771   return out;
00772 }
00773 
00774 } // namespace GLpApp
00775 
00776 //
00777 // Implementations
00778 //
00779 
00780 int GLpApp::compproduct( Epetra_SerialDenseVector & product,
00781   double *first, double *second)
00782 {
00783   for (int i=0; i<product.M(); i++) {
00784     product[i] = first[i]*second[i];
00785   }
00786   return(0);
00787 }
00788 
00789 int GLpApp::compproduct(Epetra_SerialDenseVector & product,
00790                 double *first, double *second, double *third)
00791 {
00792   for (int i=0; i<product.M(); i++) {
00793     product[i] = first[i]*second[i]*third[i];
00794   }
00795   return(0);
00796 }
00797 
00798 //#define GLPAPP_SHOW_BOUNDARY_ASSEMBLY
00799 
00800 /*  \brief Performs finite-element assembly of face mass matrices.
00801 
00802   \param  Comm      [in]  - The Epetra (MPI) communicator.
00803   \param  ipindx    [in]  - Vector of NUMIP indices of nodes that are `unique' to a subdomain
00804                             (i.e. owned by the corresponding processor).
00805   \param  ipcoords  [in]  - Matrix (NUMIP x 2) of x-y-coordinates of the vertices cooresponding
00806                             to indices ipindx: \n
00807                             ipcoords(i,0) x-coordinate of node i, \n
00808                             ipcoords(i,1) y-coordinate of node i.
00809   \param  pindx     [in]  - Vector of NUMP indices of ALL nodes in a subdomain, including
00810                             the shared nodes.
00811   \param  pcoords   [in]  - Matrix (NUMP x 2) of x-y-coordinates of the vertices cooresponding
00812                             to indices pindx: \n
00813                             pcoords(i,0) x-coordinate of node i, \n
00814                             pcoords(i,1) y-coordinate of node i.
00815   \param  t         [in]  - Matrix (ELE x 3) of indices of the vertices in a triangle: \n
00816                             t(i,j) index of the 0-based j-th vertex in triangle i, where i = 0, ..., numElements-1
00817   \param  e         [in]  - Matrix (EDGE x 3) of edges. \n
00818                             e(i,1:2) contains the indices of the endpoints of edge i, where i = 0, ..., numEdges-1 \n
00819                             e(i,3) contains the boundary marker
00820   \param  B_out     [out] - Reference-counting pointer to the Epetra_FECrsMatrix describing the FE
00821                             state/control face mass matrix.
00822   \param  R_out     [out] - Reference-counting pointer to the Epetra_FECrsMatrix describing the FE
00823                             control/control volume mass matrix.
00824   \return 0                 if successful.
00825 
00826   \par Detailed Description:
00827 
00828   -# Assembles the FE state/control mass matrix \e B, given by
00829      \f[
00830        \mathbf{B}_{jk} = b(\mu_k,\phi_j) = - \int_{\partial \Omega}  \mu_k(x) \phi_j(x) dx,
00831      \f]
00832      where \f$\{ \phi_j \}_{j = 1}^{m}\f$ is the piecewise linear nodal basis for the finite-dimensional
00833      state space, and \f$\{ \mu_j \}_{j = 1}^{n}\f$ is the piecewise linear nodal basis for the
00834      finite-dimensional control space.
00835   -# Assembles the FE control/control mass matrix \e R, given by
00836      \f[
00837        \mathbf{R}_{jk} = b(\mu_k,\mu_j) = - \int_{\partial \Omega}  \mu_k(x) \mu_j(x) dx,
00838      \f]
00839      where \f$\{ \mu_j \}_{j = 1}^{n}\f$ is the piecewise linear nodal basis for the finite-dimensional
00840      control space.
00841 */
00842 int GLpApp::assemble_bdry(
00843   const Epetra_Comm                                &Comm
00844   ,const Epetra_IntSerialDenseVector               &ipindx
00845   ,const Epetra_SerialDenseMatrix                  &ipcoords
00846   ,const Epetra_IntSerialDenseVector               &pindx
00847   ,const Epetra_SerialDenseMatrix                  &pcoords
00848   ,const Epetra_IntSerialDenseMatrix               &t
00849   ,const Epetra_IntSerialDenseMatrix               &e
00850   ,Teuchos::RefCountPtr<Epetra_FECrsMatrix>        *B_out
00851   ,Teuchos::RefCountPtr<Epetra_FECrsMatrix>        *R_out
00852   )
00853 {
00854 
00855   using Teuchos::rcp;
00856 
00857 #ifdef GLPAPP_SHOW_BOUNDARY_ASSEMBLY
00858   Teuchos::RefCountPtr<Teuchos::FancyOStream>
00859     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00860   Teuchos::OSTab tab(out);
00861   *out << "\nEntering assemble_bdry(...) ...\n";
00862 #endif
00863 
00864   int numLocElems = t.M();
00865   int numLocEdges = e.M();
00866 
00867   int indexBase = 1;
00868 
00869   //
00870   // Create a sorted (by global ID) list of boundry nodes
00871   //
00872   int * lastindx = 0;
00873   Epetra_IntSerialDenseVector BdryNode(2*numLocEdges);
00874   for (int j=0; j<numLocEdges; j++) {
00875     BdryNode[j] = e(j,0);
00876     BdryNode[j+numLocEdges] = e(j,1);
00877   }
00878   std::sort(BdryNode.Values(), BdryNode.Values()+2*numLocEdges);
00879   lastindx  = std::unique(BdryNode.Values(), BdryNode.Values()+2*numLocEdges);
00880   const int numBdryNodes = lastindx - BdryNode.Values();
00881   BdryNode.Resize(numBdryNodes); // RAB: This does not overwrite?
00882 
00883   //
00884   // Determine the boundary vertices that belong to this processor.
00885   //
00886   Epetra_IntSerialDenseVector MyBdryNode(numBdryNodes);
00887   lastindx = std::set_intersection(
00888     BdryNode.Values(), BdryNode.Values()+numBdryNodes,  // (Sorted) set A
00889     ipindx.Values(), ipindx.Values()+ipindx.M(),        // (Sorted) set B
00890     MyBdryNode.Values()                                 // Intersection
00891     );
00892   const int numMyBdryNodes = lastindx - MyBdryNode.Values();
00893   MyBdryNode.Resize(numMyBdryNodes); // RAB: This does not overwrite?
00894   
00895   //
00896   // Define the maps for the various lists
00897   //
00898   Epetra_Map standardmap(-1, ipindx.M(), const_cast<int*>(ipindx.A()), indexBase, Comm);
00899   Epetra_Map overlapmap(-1, pindx.M(), const_cast<int*>(pindx.A()), indexBase, Comm);
00900   Epetra_Map mybdryctrlmap(-1, numMyBdryNodes, const_cast<int*>(MyBdryNode.A()), indexBase, Comm);
00901   // Above, it is important to note what mybndyctrlmap represents.  It is the
00902   // a sorted list of global vertex node IDS for nodes on a boundary that are
00903   // uniquely owned by the local process.
00904 
00905 #ifdef GLPAPP_SHOW_BOUNDARY_ASSEMBLY
00906   *out << "\nstandardmap:\n";
00907   standardmap.Print(*Teuchos::OSTab(out).getOStream());
00908   *out << "\nmybdryctrlmap:\n";
00909   mybdryctrlmap.Print(*Teuchos::OSTab(out).getOStream());
00910 #endif
00911 
00912   //
00913   // Allocate matrices to fill
00914   //
00915   Teuchos::RefCountPtr<Epetra_FECrsMatrix>
00916     B = rcp(new Epetra_FECrsMatrix(Copy,standardmap,0)),
00917     R = rcp(new Epetra_FECrsMatrix(Copy,standardmap,0));
00918   // NOTE: The data map is the same as for the volume matrices. Later, when
00919   // FillComplete is called, we will fix the proper domain and range maps. 
00920   // Declare quantities needed for the call to the local assembly routine.
00921   const int numNodesPerEdge = 2;
00922   Epetra_IntSerialDenseVector epetra_nodes(numNodesPerEdge);
00923 
00924   //
00925   // Load B and R by looping through the edges
00926   //
00927 
00928   Epetra_SerialDenseMatrix Bt(2,2);
00929   int err=0;
00930 
00931   for( int i=0; i < numLocEdges; i++ ) {
00932 
00933     const int
00934       global_id_0 = e(i,0),
00935       global_id_1 = e(i,1),
00936       local_id_0  = overlapmap.LID(global_id_0), // O(log(numip)) bindary search
00937       local_id_1  = overlapmap.LID(global_id_1); // O(log(numip)) bindary search
00938 
00939     epetra_nodes(0) = global_id_0;
00940     epetra_nodes(1) = global_id_1;
00941 
00942     const double
00943       x0 = pcoords(local_id_0,0),
00944       y0 = pcoords(local_id_0,1),
00945       x1 = pcoords(local_id_1,0),
00946       y1 = pcoords(local_id_1,1);
00947     
00948     const double l = sqrt(pow(x0-x1,2) + pow(y0-y1,2));  // Length of this edge
00949     
00950     // We have an explicit formula for Bt.
00951     const double l_sixth = l * (1.0/6.0);
00952     Bt(0,0) = l_sixth * 2.0;
00953     Bt(0,1) = l_sixth * 1.0;
00954     Bt(1,0) = l_sixth * 1.0;
00955     Bt(1,1) = l_sixth * 2.0;
00956 
00957 #ifdef GLPAPP_SHOW_BOUNDARY_ASSEMBLY
00958     *out
00959       << "\nEdge global nodes = ("<<global_id_0<<","<<global_id_1<<"),"
00960       << " local nodes = ("<<local_id_0<<","<<local_id_1<<"),"
00961       << " (x0,y0)(x1,y1)=("<<x0<<","<<y0<<")("<<x1<<","<<y1<<"),"
00962       << " Bt = ["<<Bt(0,0)<<","<<Bt(0,1)<<";"<<Bt(1,0)<<","<<Bt(1,1)<<"]\n";
00963 #endif
00964 
00965     const int format = Epetra_FECrsMatrix::COLUMN_MAJOR;
00966     err = B->InsertGlobalValues(epetra_nodes,Bt,format);
00967     if (err<0) return(err);
00968     err = R->InsertGlobalValues(epetra_nodes,Bt,format);
00969     if (err<0) return(err);
00970     
00971   }
00972 
00973 /*
00974 
00975   err = B->GlobalAssemble(false);
00976   if (err<0) return(err);
00977   err = R->GlobalAssemble(false);
00978   if (err<0) return(err);
00979 
00980   err = B->FillComplete(mybdryctrlmap,standardmap);
00981   if (err<0) return(err);
00982   err = R->FillComplete(mybdryctrlmap,mybdryctrlmap);
00983   if (err<0) return(err);
00984 
00985 */
00986 
00987   err = B->GlobalAssemble(mybdryctrlmap,standardmap,true);
00988   if (err<0) return(err);
00989   err = R->GlobalAssemble(mybdryctrlmap,mybdryctrlmap,true);
00990   if (err<0) return(err);
00991 
00992   if(B_out) *B_out = B;
00993   if(R_out) *R_out = R;
00994 
00995 #ifdef GLPAPP_SHOW_BOUNDARY_ASSEMBLY
00996   *out << "\nB =\n";
00997   B->Print(*Teuchos::OSTab(out).getOStream());
00998   *out << "\nLeaving assemble_bdry(...) ...\n";
00999 #endif
01000 
01001   return(0);
01002 
01003 }
01004 
01005 /*  \brief Performs finite-element assembly of volume stiffness and mass matrices,
01006             and the right-hand side (forcing term).
01007 
01008   \param  Comm      [in]  - The Epetra (MPI) communicator.
01009   \param  ipindx    [in]  - Vector of NUMIP indices of nodes that are `unique' to a subdomain
01010                             (i.e. owned by the corresponding processor).
01011   \param  ipcoords  [in]  - Matrix (NUMIP x 2) of x-y-coordinates of the vertices cooresponding
01012                             to indices ipindx: \n
01013                             ipcoords(i,0) x-coordinate of node i, \n
01014                             ipcoords(i,1) y-coordinate of node i.
01015   \param  pindx     [in]  - Vector of NUMP indices of ALL nodes in a subdomain, including
01016                             the shared nodes.
01017   \param  pcoords   [in]  - Matrix (NUMP x 2) of x-y-coordinates of the vertices cooresponding
01018                             to indices pindx: \n
01019                             pcoords(i,0) x-coordinate of node i, \n
01020                             pcoords(i,1) y-coordinate of node i.
01021   \param  t         [in]  - Matrix (ELE x 3) of indices of the vertices in a triangle: \n
01022                             t(i,j) index of the j-th vertex in triangle i, where i = 1, ..., ELE
01023   \param  e         [in]  - Matrix (EDGE x 3) of edges. \n
01024                             e(i,1:2) contains the indices of the endpoints of edge i, where i = 1, ..., EDGE \n
01025                             e(i,3) contains the boundary marker \n
01026                             e(i,3) = 1  Dirichlet boundary conditions on edge i \n
01027                             e(i,3) = 2  Neumann boundary conditions on edge i \n
01028                             e(i,3) = 3  Robin boundary conditions on edge i
01029   \param  A         [out] - Reference-counting pointer to the Epetra_FECrsMatrix describing the FE volume
01030                             stiffness matrix for the advection-diffusion equation. Includes advection,
01031                             diffusion, and reaction terms, and modifications that account for the boundary
01032                             conditions.
01033   \param  M         [out] - Reference-counting pointer to the Epetra_FECrsMatrix describing the FE volume
01034                             mass matrix.
01035   \param  b         [out] - Reference-counting pointer to the Epetra_FEVector describing the discretized
01036                             forcing term in the advection-diffusion equation. Includes modifications that
01037                             account for the boundary conditions.
01038   \return 0                 if successful.
01039 
01040   \par Detailed Description:
01041 
01042   -# Assembles the FE volume stiffness matrix \e A and the right-hand side \e b for the
01043   solution of an advection-diffusion equation using piecewise linear finite elements.
01044   The advection-diffusion equation is 
01045   \f{align*}
01046        - \nabla \cdot \left( K(x) \nabla y(x) \right) + \mathbf{c}(x) \cdot \nabla y(x) + r(x) y(x)
01047        &= f(x), & x &\in \Omega,\;  \\
01048        y(x) &= d(x),    & x &\in {\partial \Omega}_d, \\
01049        K(x) \frac{\partial}{\partial \mathbf{n}}  y(x) &= g(x), & x &\in {\partial \Omega}_n, \\
01050        \sigma_0 K(x) \frac{\partial}{\partial \mathbf{n}} y(x)
01051        + \sigma_1 y(x) &= g(x), & x &\in {\partial \Omega}_r,
01052   \f}
01053   where \f$ K \f$ represents scalar diffusion, \f$ \mathbf{c} \f$ is the advection vector field,
01054   \f$ r \f$ is the reaction term, \f$ d \f$ and  \f$ g \f$ are given functions, \f$sigma_0\f$ and
01055   \f$ sigma_1 \f$ are given quantities, \f$ {\partial \Omega}_d \f$ is the Dirichlet boundary,
01056   \f$ {\partial \Omega}_n \f$ is the Neumann boundary, and \f$ {\partial \Omega}_r \f$ is the Robin
01057   boundary. The quantities \f$ K \f$, \f$ \mathbf{c} \f$, \f$ r \f$, \f$ d \f$, and \f$ g \f$ are
01058   assumed to be piecewise linear. Currently, they are to be hard-coded inside this function.
01059   -# Assembles the FE volume mass matrix \e M.
01060 */
01061 int GLpApp::assemble(const Epetra_Comm & Comm,
01062              const Epetra_IntSerialDenseVector & ipindx,
01063              const Epetra_SerialDenseMatrix & ipcoords,
01064              const Epetra_IntSerialDenseVector & pindx,
01065              const Epetra_SerialDenseMatrix & pcoords,
01066              const Epetra_IntSerialDenseMatrix & t,
01067              const Epetra_IntSerialDenseMatrix & e,
01068              RefCountPtr<Epetra_FECrsMatrix> & A,
01069              RefCountPtr<Epetra_FECrsMatrix> & M,
01070              RefCountPtr<Epetra_FEVector> & b)
01071 {
01072 
01073   int myPID = Comm.MyPID();
01074   int numProcs = Comm.NumProc();
01075   Usr_Par usr_par;
01076 
01077   int numLocElems = t.M();
01078   int numNodesPerElem = 3;
01079 
01080   int indexBase = 1;
01081 
01082   Epetra_Map standardmap(-1, ipindx.M(), (int*)ipindx.A(), indexBase, Comm);
01083   Epetra_Map overlapmap(-1, pindx.M(), (int*)pindx.A(), indexBase, Comm);
01084 
01085   int* nodes = new int[numNodesPerElem];
01086   int i=0, j=0, err=0;
01087 
01088   A = rcp(new Epetra_FECrsMatrix(Copy, standardmap, 0));
01089   M = rcp(new Epetra_FECrsMatrix(Copy, standardmap, 0));
01090   b = rcp(new Epetra_FEVector(standardmap));
01091 
01092   // Declare quantities needed for the call to the local assembly routine.
01093   int format = Epetra_FECrsMatrix::COLUMN_MAJOR;
01094   Epetra_IntSerialDenseVector epetra_nodes(View, nodes, numNodesPerElem);
01095 
01096 
01097   /* ************************  First, we build A and b.  ************************ */
01098   Epetra_SerialDenseMatrix At;
01099   Epetra_SerialDenseVector bt;
01100   Epetra_SerialDenseMatrix vertices(numNodesPerElem, pcoords.N());
01101   
01102   Epetra_SerialDenseVector k(numNodesPerElem);
01103   for (i=0; i< numNodesPerElem; i++) k(i)=1.0;
01104   Epetra_SerialDenseMatrix c(numNodesPerElem,2);
01105   for (i=0; i< numNodesPerElem; i++) { c(i,0)=0.0; c(i,1)=0.0; }
01106   Epetra_SerialDenseVector r(numNodesPerElem);
01107   for (i=0; i< numNodesPerElem; i++) r(i)=0.0;
01108   Epetra_SerialDenseVector f(numNodesPerElem);
01109   for (i=0; i< numNodesPerElem; i++) f(i)=0.0;
01110   Epetra_SerialDenseVector g(2);
01111   g(0)=0.0; g(1)=0.0;
01112   Epetra_SerialDenseVector sigma(2);
01113   sigma(0)=0.0; sigma(1)=0.0;
01114   for(i=0; i<numLocElems; i++) {
01115     nodes[0] = t(i,0); nodes[1] = t(i,1); nodes[2] = t(i,2);
01116     for (j=0; j<numNodesPerElem; j++) {
01117       // get vertex
01118       vertices(j,0) = pcoords(overlapmap.LID(nodes[j]), 0);
01119       vertices(j,1) = pcoords(overlapmap.LID(nodes[j]), 1);
01120       // set rhs function
01121       f(j) = cos(GLp_pi*vertices(j,0))*cos(GLp_pi*vertices(j,1)) * 
01122                (2*GLp_pi*GLp_pi + pow(cos(GLp_pi*vertices(j,0)),2)*pow(cos(GLp_pi*vertices(j,1)),2) - 1.0);
01123     }
01124     lassembly(vertices, k, c, r, f, usr_par, At, bt);
01125     err = A->InsertGlobalValues(epetra_nodes, At, format);
01126     if (err<0) return(err);
01127     err = b->SumIntoGlobalValues(epetra_nodes, bt);
01128     if (err<0) return(err);
01129   }
01130 
01131   // MAKE ADJUSTMENTS TO A and b TO ACCOUNT FOR BOUNDARY CONDITIONS.
01132 
01133   // Get Neumann boundary edges.
01134   Epetra_IntSerialDenseMatrix neumann(e.M(),2);
01135   j = 0;
01136   for (i=0; i<e.M(); i++) {
01137     if (e(i,2)==2) {
01138       neumann(j,0) = e(i,0);  neumann(j,1) = e(i,1);
01139       j++;
01140     }
01141   }
01142   neumann.Reshape(j,2);
01143   // Adjust for Neumann BC's.
01144   if (neumann.M() != 0) {
01145     Epetra_BLAS blas;
01146     Epetra_SerialDenseMatrix quadnodes;
01147     Epetra_SerialDenseVector quadweights;
01148     Epetra_SerialDenseMatrix N;
01149     Epetra_SerialDenseMatrix NN;
01150     Epetra_SerialDenseVector product(2);
01151 
01152     // Get quadrature weights and points.
01153     quadrature(1, 2, quadnodes, quadweights);
01154     
01155     // Evaluate nodal basis functions at quadrature points
01156     // N(i,j) value of basis function j at quadrature node i
01157     N.Shape(quadnodes.M(),2);
01158     for (i=0; i<quadnodes.M(); i++) {
01159       N(i,0) = 1.0 - quadnodes(i,0);
01160       N(i,1) = quadnodes(i,0);
01161     }
01162 
01163     // Evaluate integrals of 4 products N(i)*N(j).
01164     NN.Shape(2,2);
01165     for (i=0; i<2; i++) {
01166       for (j=0; j<2; j++) {
01167         compproduct(product, N[i], N[j]);
01168         NN(i,j) = blas.DOT(quadweights.M(), quadweights.A(), product.A());
01169       }
01170     }
01171 
01172     Epetra_IntSerialDenseVector neumann_nodes(2);
01173     Epetra_SerialDenseVector neumann_add(2);
01174     double l;
01175     for (i=0; i<neumann.M(); i++) {
01176       neumann_nodes(0) = neumann(i,0); neumann_nodes(1) = neumann(i,1);
01177       neumann_add(0) = pcoords(overlapmap.LID(neumann_nodes(0)),0)
01178                       -pcoords(overlapmap.LID(neumann_nodes(1)),0);
01179       neumann_add(1) = pcoords(overlapmap.LID(neumann_nodes(0)),1)
01180                       -pcoords(overlapmap.LID(neumann_nodes(1)),1);
01181       l = blas.NRM2(neumann_add.M(), neumann_add.A());
01182       neumann_add.Multiply('N', 'N', 1.0, NN, g, 0.0);
01183       neumann_add.Scale(l);
01184       err = b->SumIntoGlobalValues(neumann_nodes, neumann_add); 
01185       if (err<0) return(err);
01186     }
01187   }
01188 
01189   // Get Robin boundary edges.
01190   Epetra_IntSerialDenseMatrix robin(e.M(),2);
01191   j = 0;
01192   for (i=0; i<e.M(); i++) {
01193     if (e(i,2)==3) {
01194       robin(j,0) = e(i,0);  robin(j,1) = e(i,1);
01195       j++;
01196     }
01197   }
01198   robin.Reshape(j,2);
01199   // Adjust for Robin BC's.
01200   if (robin.M() != 0) {
01201     Epetra_BLAS blas;
01202     Epetra_SerialDenseMatrix quadnodes;
01203     Epetra_SerialDenseVector quadweights;
01204     Epetra_SerialDenseMatrix N;
01205     Epetra_SerialDenseMatrix NN;
01206     Epetra_SerialDenseMatrix * NNN;
01207     Epetra_SerialDenseVector product(2);
01208 
01209     // Get quadrature weights and points.
01210     quadrature(1, 2, quadnodes, quadweights);
01211     
01212     // Evaluate nodal basis functions at quadrature points
01213     // N(i,j) value of basis function j at quadrature node i
01214     N.Shape(quadnodes.M(),2);
01215     for (i=0; i<quadnodes.M(); i++) {
01216       N(i,0) = 1.0 - quadnodes(i,0);
01217       N(i,1) = quadnodes(i,0);
01218     }
01219 
01220     // Evaluate integrals of 4 products N(i)*N(j).
01221     NN.Shape(2,2);
01222     for (i=0; i<2; i++) {
01223       for (j=0; j<2; j++) {
01224         compproduct(product, N[i], N[j]);
01225         NN(i,j) = blas.DOT(quadweights.M(), quadweights.A(), product.A());
01226       }
01227     }
01228 
01229     // Evaluate integrals of 8 products N(i)*N(j)*N(k).
01230     NNN = new Epetra_SerialDenseMatrix[2];
01231     NNN[0].Shape(2,2); NNN[1].Shape(2,2);
01232     for (i=0; i<2; i++) {
01233       for (j=0; j<2; j++) {
01234         for (int k=0; k<2; k++) {
01235           compproduct(product, N[i], N[j], N[k]);
01236           NNN[k](i,j) = blas.DOT(quadweights.M(), quadweights.A(),
01237                                  product.A());
01238         }
01239       }
01240     }
01241 
01242     Epetra_IntSerialDenseVector robin_nodes(2);
01243     Epetra_SerialDenseVector robin_b_add(2);
01244     Epetra_SerialDenseMatrix robin_A_add(2,2);
01245     double l;
01246     for (i=0; i<robin.M(); i++) {
01247       robin_nodes(0) = robin(i,0); robin_nodes(1) = robin(i,1);
01248       
01249       robin_b_add(0) = pcoords(overlapmap.LID(robin_nodes(0)),0)
01250                       -pcoords(overlapmap.LID(robin_nodes(1)),0);
01251       robin_b_add(1) = pcoords(overlapmap.LID(robin_nodes(0)),1)
01252                       -pcoords(overlapmap.LID(robin_nodes(1)),1);
01253       l = blas.NRM2(robin_b_add.M(), robin_b_add.A());
01254       robin_b_add.Multiply('N', 'N', 1.0, NN, g, 0.0);
01255       robin_b_add.Scale(l);
01256       err = b->SumIntoGlobalValues(robin_nodes, robin_b_add); 
01257       if (err<0) return(err);
01258 
01259       NNN[0].Scale(sigma(0)); NNN[1].Scale(sigma(1));
01260       robin_A_add += NNN[0]; robin_A_add += NNN[1];
01261       robin_A_add.Scale(l);
01262       err = A->InsertGlobalValues(robin_nodes, robin_A_add, format);
01263       if (err<0) return(err);
01264     }
01265 
01266     delete [] NNN;
01267   }
01268 
01269   // Get Dirichlet boundary edges.
01270   Epetra_IntSerialDenseMatrix dirichlet(e.M(),2);
01271   j = 0;
01272   for (i=0; i<e.M(); i++) {
01273     if (e(i,2)==1) {
01274       dirichlet(j,0) = e(i,0);  dirichlet(j,1) = e(i,1);
01275       j++;
01276     }
01277   }
01278   dirichlet.Reshape(j,2);
01279   // DIRICHLET NOT DONE! DO THIS LATER!!!!
01280 
01281   /* ************************  Done building A and b.  ************************ */
01282 
01283 
01284 
01285   /* ************************  Second, we build M.  ************************ */
01286 
01287   Epetra_SerialDenseMatrix Mt;
01288 
01289   for (i=0; i< numNodesPerElem; i++) k(i)=0.0;
01290   for (i=0; i< numNodesPerElem; i++) { c(i,0)=0.0; c(i,1)=0.0; }
01291   for (i=0; i< numNodesPerElem; i++) r(i)=1.0;
01292   for (i=0; i< numNodesPerElem; i++) f(i)=0.0;
01293   g(0)=0.0; g(1)=0.0;
01294   sigma(0)=0.0; sigma(1)=0.0;
01295   for(i=0; i<numLocElems; i++) {
01296     nodes[0] = t(i,0); nodes[1] = t(i,1); nodes[2] = t(i,2);
01297     for (j=0; j<numNodesPerElem; j++) {
01298       vertices(j,0) = pcoords(overlapmap.LID(nodes[j]), 0);
01299       vertices(j,1) = pcoords(overlapmap.LID(nodes[j]), 1);
01300     }
01301     lassembly(vertices, k, c, r, f, usr_par, Mt, bt);
01302     err = M->InsertGlobalValues(epetra_nodes, Mt, format);
01303     if (err<0) return(err);
01304   }
01305 
01306   /* ************************  Done building M.  ************************ */
01307 
01308 
01309 
01310   // Call global assemble and FillComplete at the same time.
01311 
01312   err = A->GlobalAssemble();
01313   if (err<0) return(err);
01314 
01315   err = b->GlobalAssemble();
01316   if (err<0) return(err);
01317 
01318   err = M->GlobalAssemble();
01319   if (err<0) return(err);
01320 
01321   delete [] nodes;
01322 
01323   return(0);
01324 }
01325 
01326 /* \brief Computes local stiffness matrix and local RHS vector for simplex
01327            (triangles in two dimensions).
01328                             
01329   \param  vertices   [in]  - Matrix containing the global coordinates of the current simplex.
01330   \param  k          [in]  - Vector containing the value of the diffusion \f$k\f$ at each vertex
01331                              (\f$k\f$ is assumed to be piecewise linear), where \f$k\f$ is the coeff in the
01332                              term \f$ \nabla \cdot (k \nabla(u)) \f$ in the advection-diffusion equation.
01333   \param  c          [in]  - Matrix containing the value of the advection field \f$ \mathbf{c} \f$ at each
01334                              vertex (\f$ \mathbf{c} \f$ is assumed to be piecewise linear), where
01335                              \f$ \mathbf{c} \f$ is the 2-vector of coeffs in the term
01336                              \f$ \mathbf{c}(x) \cdot \nabla y(x) \f$ in the advection-diffusion equation.
01337   \param  r          [in]  - Vector containing the value of \f$ r \f$ at each vertex (\f$ r \f$ is assumed
01338                              to be piecewise linear), where \f$ r \f$ is the coefficient in the term
01339                              \f$ ru \f$ in the advection-diffusion equation.
01340   \param  f          [in]  - Vector containing the value of \f$ f \f$ at each vertex (\f$ f \f$ is assumed to be
01341                              piecewise linear), where \f$ f \f$ is the right hand side in the adv-diff eq
01342   \param  usr_par    [in]  - class containing: \n
01343                               - S1, S2, S3 (3x3) the integrals of various combinations of partials
01344                                 of local basis functions
01345                               - N (1x3) integrals of local basis functions
01346                               - NNN[3] (3x3), integrals of products of local basis functions N(i)*N(j)*N(k)
01347                               - etc.
01348   \param  At         [out] - stiffness matrix contribution for the simplex
01349   \param  bt         [out] - right-hand side contribution for the simplex
01350 
01351   \return 0                 if successful.
01352 
01353   \par Detailed Description:
01354 
01355      Computes the local (per simplex) contributions to the FE volume stiffness matrix \e A and the
01356      right-hand side \e b for the advection-diffusion equation
01357     \f{align*}
01358        - \nabla \cdot \left( K(x) \nabla y(x) \right) + \mathbf{c}(x) \cdot \nabla y(x) + r(x) y(x)
01359        &= f(x), & x &\in \Omega,\;  \\
01360        y(x) &= d(x),    & x &\in {\partial \Omega}_d, \\
01361        K(x) \frac{\partial}{\partial \mathbf{n}}  y(x) &= g(x), & x &\in {\partial \Omega}_n, \\
01362        \sigma_0 K(x) \frac{\partial}{\partial \mathbf{n}} y(x)
01363        + \sigma_1 y(x) &= g(x), & x &\in {\partial \Omega}_r.
01364     \f}
01365 */
01366 int GLpApp::lassembly(const Epetra_SerialDenseMatrix & vertices,
01367               const Epetra_SerialDenseVector & k,
01368               const Epetra_SerialDenseMatrix & c,
01369               const Epetra_SerialDenseVector & r,
01370               const Epetra_SerialDenseVector & f,
01371               const Usr_Par & usr_par,
01372               Epetra_SerialDenseMatrix & At,
01373               Epetra_SerialDenseVector & bt)
01374 {
01375   // Note that the constructors below initialize all entries to 0.
01376   Epetra_SerialDenseMatrix B(2,2);
01377   Epetra_SerialDenseMatrix b(2,2);
01378   Epetra_SerialDenseMatrix BtB(2,2);  
01379   Epetra_SerialDenseMatrix C(2,2);  
01380   Epetra_SerialDenseMatrix M1(3,3);
01381   Epetra_SerialDenseMatrix M2(3,3);
01382   Epetra_SerialDenseMatrix M3(3,3);
01383   Epetra_SerialDenseMatrix tmp(3,3);
01384   double dB, adB;
01385   At.Shape(3,3);
01386   bt.Size(3);
01387 
01388   // Construct affine transformation matrix.
01389   for(int i=0; i<2; i++) {
01390     B(i,0) = vertices(1,i)-vertices(0,i);
01391     B(i,1) = vertices(2,i)-vertices(0,i);
01392   }
01393   dB  = determinant(B);
01394   adB = abs(dB);
01395 
01396   // Compute matrix C = inv(B'*B).
01397   BtB.Multiply('T', 'N', 1.0, B, B, 0.0);
01398   inverse(BtB, C);
01399 
01400   inverse(B, b); b.Scale(dB);
01401   
01402   // Compute the part corresponding to div(K*grad(u)).
01403   tmp = usr_par.S1; tmp.Scale(C(0,0));
01404   M1 += tmp;
01405   tmp = usr_par.S2; tmp.Scale(C(0,1));
01406   M1 += tmp;
01407   tmp = usr_par.S3; tmp.Scale(C(1,1));
01408   M1 += tmp;
01409   M1.Scale( (k(0)*usr_par.Nw(0) + k(1)*usr_par.Nw(1) +
01410              k(2)*usr_par.Nw(2)) * adB );
01411 
01412   // Compute the part corresponding to c'*grad(u).
01413   tmp = usr_par.NdNdx1Nw[0]; tmp.Scale(b(0,0)*c(0,0)+b(0,1)*c(0,1));
01414   M2 += tmp;
01415   tmp = usr_par.NdNdx1Nw[1]; tmp.Scale(b(0,0)*c(1,0)+b(0,1)*c(1,1));
01416   M2 += tmp;
01417   tmp = usr_par.NdNdx1Nw[2]; tmp.Scale(b(0,0)*c(2,0)+b(0,1)*c(2,1));
01418   M2 += tmp;
01419   tmp = usr_par.NdNdx2Nw[0]; tmp.Scale(b(1,0)*c(0,0)+b(1,1)*c(0,1));
01420   M2 += tmp;
01421   tmp = usr_par.NdNdx2Nw[1]; tmp.Scale(b(1,0)*c(1,0)+b(1,1)*c(1,1));
01422   M2 += tmp;
01423   tmp = usr_par.NdNdx2Nw[2]; tmp.Scale(b(1,0)*c(2,0)+b(1,1)*c(2,1));
01424   M2 += tmp;
01425   M2.Scale(adB/dB);
01426 
01427   // Compute the part corresponding to r*u.
01428   tmp = usr_par.NNNw[0]; tmp.Scale(r(0));
01429   M3 += tmp;
01430   tmp = usr_par.NNNw[1]; tmp.Scale(r(1));
01431   M3 += tmp;
01432   tmp = usr_par.NNNw[2]; tmp.Scale(r(2));
01433   M3 += tmp;
01434   M3.Scale(adB);
01435 
01436   At += M1;
01437   At += M2;
01438   At += M3;
01439 
01440   bt.Multiply('T', 'N', adB, usr_par.NNw, f, 0.0);  
01441   
01442   return(0);
01443 }
01444 
01445 /*  \brief Computes the inverse of a dense matrix.
01446 
01447   \param  mat  [in]  - the matrix that is to be inverted
01448   \param  inv  [in]  - the inverse
01449 
01450   \return 0            if successful
01451 */
01452 int GLpApp::inverse(const Epetra_SerialDenseMatrix & mat,
01453             Epetra_SerialDenseMatrix & inv)
01454 {
01455   Epetra_LAPACK lapack;
01456   int dim = mat.M();
01457   int info;
01458   Epetra_IntSerialDenseVector ipiv(dim);
01459   Epetra_SerialDenseVector work(2*dim);
01460 
01461   inv.Shape(dim,dim);
01462   inv = mat;
01463 
01464   lapack.GETRF(dim, dim, inv.A(), dim, ipiv.A(), &info);
01465   lapack.GETRI(dim, inv.A(), dim, ipiv.A(), work.A(), &dim, &info);
01466   
01467   return(0);
01468 }
01469 
01470 /*  \brief Computes the determinant of a dense matrix.
01471 
01472   \param  mat  [in]  - the matrix
01473 
01474   \return the determinant
01475 */
01476 double GLpApp::determinant(const Epetra_SerialDenseMatrix & mat)
01477 {
01478   //Teuchos::LAPACK<int,double> lapack;
01479   Epetra_LAPACK lapack;
01480   double det;
01481   int swaps; 
01482   int dim = mat.M();
01483   int info;
01484   Epetra_IntSerialDenseVector ipiv(dim);
01485  
01486   Epetra_SerialDenseMatrix mymat(mat);
01487   lapack.GETRF(dim, dim, mymat.A(), dim, ipiv.A(), &info);
01488 
01489   det = 1.0;
01490   for (int i=0; i<dim; i++) {
01491     det *= mymat(i,i);
01492   }
01493 
01494   swaps = 0;
01495   for (int i=0; i<dim; i++) {
01496     if ((ipiv[i]-1) != i)
01497       swaps++;
01498   }
01499 
01500   det *= pow((double)(-1.0),swaps);
01501 
01502   return(det);
01503 }
01504 
01505 int GLpApp::meshreader(const Epetra_Comm & Comm,
01506                Epetra_IntSerialDenseVector & ipindx,
01507                Epetra_SerialDenseMatrix & ipcoords,
01508                Epetra_IntSerialDenseVector & pindx,
01509                Epetra_SerialDenseMatrix & pcoords,
01510                Epetra_IntSerialDenseMatrix & t,
01511                Epetra_IntSerialDenseMatrix & e,
01512                const char geomFileBase[],
01513                const bool trace,
01514                const bool dumpAll
01515                )
01516 {
01517   int MyPID = Comm.MyPID();
01518 
01519   const int FileNameSize = 120;
01520   char FileName[FileNameSize];
01521   TEST_FOR_EXCEPT(static_cast<int>(std::strlen(geomFileBase) + 5) > FileNameSize);
01522   sprintf(FileName, "%s.%03d", geomFileBase, MyPID);
01523 
01524   if(1) {
01525     std::ifstream file_in(FileName);
01526     TEST_FOR_EXCEPTION(
01527       file_in.eof(), std::logic_error
01528       ,"Error, the file \""<<FileName<<"\" could not be opened for input!"
01529       );
01530   }
01531 
01532   FILE* fid;
01533   fid = fopen(FileName, "r");
01534 
01535   if(trace) printf("\nReading node info from %s ...\n", FileName);
01536   int numip = 0, numcp = 0; // # owned nodes, # shared nodes
01537   fscanf(fid, "%d %d", &numip, &numcp);
01538   const int nump = numip + numcp; // # total nodes
01539   if(trace) printf( "\nnumip = %d, numcp = %d, nump = %d\n", numip, numcp, nump );
01540   ipindx.Size(numip);
01541   ipcoords.Shape(numip, 2);
01542   pindx.Size(nump);
01543   pcoords.Shape(nump, 2);
01544   for (int i=0; i<numip; i++) {
01545     fscanf(fid,"%d %lf %lf %*d",&ipindx(i),&ipcoords(i,0),&ipcoords(i,1));
01546     if(trace&&dumpAll) printf("%d %lf %lf\n",ipindx(i),ipcoords(i,0),ipcoords(i,1));
01547     pindx(i) = ipindx(i);
01548     pcoords(i,0) = ipcoords(i,0); pcoords(i,1) = ipcoords(i,1);
01549   }
01550   for (int i=numip; i<nump; i++) {
01551     fscanf(fid,"%d %lf %lf %*d",&pindx(i),&pcoords(i,0),&pcoords(i,1));
01552   }
01553 
01554   fscanf(fid,"%*[^\n]");   // Skip to the End of the Line
01555   fscanf(fid,"%*1[\n]");   // Skip One Newline
01556 
01557   fscanf(fid,"%*[^\n]");   // Skip to the End of the Line
01558   fscanf(fid,"%*1[\n]");   // Skip One Newline
01559 
01560   for (int i=0; i<nump; i++) {
01561     fscanf(fid,"%*[^\n]"); // Skip to the End of the Line
01562     fscanf(fid,"%*1[\n]"); // Skip One Newline
01563   }
01564 
01565   if(trace) printf("\nReading element info from %s ...\n", FileName);
01566   int numelems = 0; // # elements on this processor
01567   fscanf(fid, "%d", &numelems);
01568   if(trace) printf( "\nnumelems = %d\n", numelems );
01569   t.Shape(numelems,3);
01570   for (int i=0; i<numelems; i++) {
01571     fscanf(fid, "%d %d %d", &t(i,0), &t(i,1), &t(i,2));
01572     if(trace&&dumpAll) printf("%d %d %d\n", t(i,0), t(i,1), t(i,2));
01573   }
01574 
01575   if(trace) printf("\nReading boundry edge info from %s ...\n", FileName);
01576   int numedges = 0; // # boundry edges on this processor
01577   fscanf(fid,"%d",&numedges);
01578   if(trace) printf( "\nnumedges = %d\n", numedges );
01579   e.Shape(numedges,3);
01580   for (int i=0; i<numedges; i++) {
01581     fscanf(fid, "%d %d %d", &e(i,0), &e(i,1), &e(i,2));
01582     if(trace&&dumpAll) printf("%d %d %d\n", e(i,0), e(i,1), e(i,2));
01583   }
01584 
01585   fclose(fid);
01586   if(trace) printf("Done reading mesh.\n");
01587 
01588   return(0);
01589 
01590 }
01591 
01592 /*  \brief Performs finite-element assembly of the Hessian of the nonlinear form times an arbitrary vector.
01593 
01594   \param  Comm      [in]  - The Epetra (MPI) communicator.
01595   \param  ipindx    [in]  - Vector of NUMIP indices of nodes that are `unique' to a subdomain
01596                             (i.e. owned by the corresponding processor).
01597   \param  ipcoords  [in]  - Matrix (NUMIP x 2) of x-y-coordinates of the vertices cooresponding
01598                             to indices ipindx: \n
01599                             ipcoords(i,0) x-coordinate of node i, \n
01600                             ipcoords(i,1) y-coordinate of node i.
01601   \param  pindx     [in]  - Vector of NUMP indices of ALL nodes in a subdomain, including
01602                             the shared nodes.
01603   \param  pcoords   [in]  - Matrix (NUMP x 2) of x-y-coordinates of the vertices cooresponding
01604                             to indices pindx: \n
01605                             pcoords(i,0) x-coordinate of node i, \n
01606                             pcoords(i,1) y-coordinate of node i.
01607   \param  t         [in]  - Matrix (ELE x 3) of indices of the vertices in a triangle: \n
01608                             t(i,j) index of the j-th vertex in triangle i, where i = 1, ..., ELE
01609   \param  y         [in]  - Reference-counting pointer to the Epetra_MultiVector at which the nonlinear
01610                             term is evaluated.
01611   \param  s         [in]  - Reference-counting pointer to the Epetra_MultiVector that is multiplied by the
01612                             Hessian of the nonlinear form.
01613   \param  lambda    [in]  - Reference-counting pointer to the Epetra_MultiVector of Lagrange Multipliers.
01614   \param  hessvec   [out] - Reference-counting pointer to the Epetra_FECrsMatrix containing Hessian of
01615                             the nonlinear form times vector product.
01616   \return 0                 if successful.
01617 
01618   \par Detailed Description:
01619 
01620   Assembles the nonlinear term \e hessvec, represented by
01621      \f[
01622        \{N''(y,\lambda)s\}_{j} = \langle g''(y_h) \lambda s,\phi_j \rangle
01623                                = \int_{\Omega} g''(y_h(x)) \lambda(x) s(x) \phi_j(x) dx,
01624      \f]
01625   where \f$ g''(y_h) \f$ is given in the local function \e g2pfctn, and \f$\{ \phi_j \}_{j = 1}^{m}\f$ is the
01626   piecewise linear nodal basis for the state space.
01627 */
01628 int GLpApp::nonlinhessvec(const Epetra_Comm & Comm,
01629                   const Epetra_IntSerialDenseVector & ipindx,
01630                   const Epetra_SerialDenseMatrix & ipcoords,
01631                   const Epetra_IntSerialDenseVector & pindx,
01632                   const Epetra_SerialDenseMatrix & pcoords,
01633                   const Epetra_IntSerialDenseMatrix & t,
01634                   const Teuchos::RefCountPtr<const Epetra_MultiVector> & y,
01635                   const Teuchos::RefCountPtr<const Epetra_MultiVector> & s,
01636                   const Teuchos::RefCountPtr<const Epetra_MultiVector> & lambda,
01637                   Teuchos::RefCountPtr<Epetra_FEVector> & hessvec)
01638 {
01639 
01640   int myPID = Comm.MyPID();
01641   int numProcs = Comm.NumProc();
01642 
01643   int numLocNodes     = pindx.M();
01644   int numMyLocNodes   = ipindx.M();
01645   int numLocElems     = t.M();
01646   int numNodesPerElem = 3;
01647 
01648   int indexBase = 1;
01649 
01650   Epetra_Map standardmap(-1, numMyLocNodes, (int*)ipindx.A(), indexBase, Comm);
01651   Epetra_Map overlapmap(-1, numLocNodes, (int*)pindx.A(), indexBase, Comm);
01652 
01653   hessvec = rcp(new Epetra_FEVector(standardmap));
01654 
01655   int* nodes = new int[numNodesPerElem];
01656   int i=0, j=0, err=0;
01657   
01658   // get quadrature nodes and weights
01659   Epetra_SerialDenseMatrix Nodes;
01660   Epetra_SerialDenseVector Weights;
01661   quadrature(2,3,Nodes,Weights);
01662   int numQuadPts = Nodes.M();
01663 
01664   // Evaluate nodal basis functions and their derivatives at quadrature points
01665   // N(i,j) = value of the j-th basis function at quadrature node i.
01666   Epetra_SerialDenseMatrix N;
01667   N.Shape(numQuadPts,3);
01668   for (int i=0; i<numQuadPts; i++) {
01669     N(i,0) = 1.0 - Nodes(i,0) - Nodes(i,1);
01670     N(i,1) = Nodes(i,0);
01671     N(i,2) = Nodes(i,1);
01672   }
01673 
01674   // Declare quantities needed for the call to the local assembly routine.
01675   Epetra_IntSerialDenseVector epetra_nodes(View, nodes, numNodesPerElem);
01676   Epetra_SerialDenseMatrix vertices(numNodesPerElem, pcoords.N());
01677 
01678   Epetra_SerialDenseVector ly;        // local entries of y
01679   Epetra_SerialDenseVector Nly;       // N*ly
01680   Epetra_SerialDenseVector ls;        // local entries of s
01681   Epetra_SerialDenseVector Nls;       // N*ls
01682   Epetra_SerialDenseVector llambda;   // local entries of lambda
01683   Epetra_SerialDenseVector Nllambda;  // N*llambda
01684   Epetra_SerialDenseVector lgfctn;    // gfctn(Nly)
01685   Epetra_SerialDenseVector lgfctnNi;  // lgfctn.*N(:,i)
01686   Epetra_SerialDenseVector lgfctnNls; // lgfctn.*N(:,i).*llambda.*ls
01687   Epetra_SerialDenseVector lhessvec;  // local contribution
01688   // Size and init to zero.
01689   ly.Size(numNodesPerElem);
01690   Nly.Size(numQuadPts);
01691   ls.Size(numNodesPerElem);
01692   Nls.Size(numQuadPts);
01693   llambda.Size(numNodesPerElem);
01694   Nllambda.Size(numQuadPts);
01695   lgfctn.Size(numQuadPts);
01696   lgfctnNi.Size(numQuadPts);
01697   lgfctnNls.Size(numQuadPts);
01698   lhessvec.Size(numNodesPerElem);
01699   
01700   Epetra_SerialDenseMatrix B(2,2);
01701   double adB;
01702   
01703   for(i=0; i<numLocElems; i++) {
01704 
01705     nodes[0] = t(i,0); nodes[1] = t(i,1); nodes[2] = t(i,2);
01706     for (j=0; j<numNodesPerElem; j++) {
01707       vertices(j,0) = pcoords(overlapmap.LID(nodes[j]), 0);
01708       vertices(j,1) = pcoords(overlapmap.LID(nodes[j]), 1);
01709     }
01710 
01711     // Construct affine transformation matrix.
01712     for(int i=0; i<2; i++) {
01713       B(i,0) = vertices(1,i)-vertices(0,i);
01714       B(i,1) = vertices(2,i)-vertices(0,i);
01715     }
01716     adB  = abs(determinant(B));
01717 
01718     // Construct local (to each processor) element view of y, s, l.
01719     for (j=0; j<numNodesPerElem; j++) {
01720       ly(j) = (*((*y)(0)))[overlapmap.LID(nodes[j])];
01721       ls(j) = (*((*s)(0)))[overlapmap.LID(nodes[j])];
01722       llambda(j) = (*((*lambda)(0)))[overlapmap.LID(nodes[j])];
01723     }
01724 
01725     Nly.Multiply('N', 'N', 1.0, N, ly, 0.0);            //N*ly
01726     Nls.Multiply('N', 'N', 1.0, N, ls, 0.0);            //N*ls
01727     Nllambda.Multiply('N', 'N', 1.0, N, llambda, 0.0);  //N*llambda
01728     g2pfctn(Nly, lgfctn);
01729     
01730     for (int i=0; i<numNodesPerElem; i++) {
01731       compproduct(lgfctnNi, lgfctn.A(), N[i]);
01732       compproduct(lgfctnNls, lgfctnNi.A(), Nls.A(), Nllambda.A());  // g2pfctn(Nly).*Nls.*Nllambda.*N(:,i)
01733       lhessvec(i) = adB*lgfctnNls.Dot(Weights);
01734     }
01735 
01736     err = hessvec->SumIntoGlobalValues(epetra_nodes, lhessvec);
01737     if (err<0) return(err);
01738   }
01739 
01740   // Call global assemble.
01741 
01742   err = hessvec->GlobalAssemble();
01743   if (err<0) return(err);
01744 
01745   delete [] nodes;
01746 
01747   return(0);
01748 }
01749 
01750 /*  \brief Componentwise evaluation of the second derivative of the nonlinear reaction term.
01751   \param  v   [in]  - Vector at which the second derivative is evaluated.
01752   \param  gv  [out] - Vector value.
01753 */
01754 void GLpApp::g2pfctn(const Epetra_SerialDenseVector & v, Epetra_SerialDenseVector & gv) {
01755   for (int i=0; i<v.M(); i++) {
01756     gv(i) = 6.0*v(i);
01757   }  
01758 }
01759 
01760 /*  \brief Performs finite-element assembly of the Jacobian of the nonlinear form.
01761 
01762   \param  Comm      [in]  - The Epetra (MPI) communicator.
01763   \param  ipindx    [in]  - Vector of NUMIP indices of nodes that are `unique' to a subdomain
01764                             (i.e. owned by the corresponding processor).
01765   \param  ipcoords  [in]  - Matrix (NUMIP x 2) of x-y-coordinates of the vertices cooresponding
01766                             to indices ipindx: \n
01767                             ipcoords(i,0) x-coordinate of node i, \n
01768                             ipcoords(i,1) y-coordinate of node i.
01769   \param  pindx     [in]  - Vector of NUMP indices of ALL nodes in a subdomain, including
01770                             the shared nodes.
01771   \param  pcoords   [in]  - Matrix (NUMP x 2) of x-y-coordinates of the vertices cooresponding
01772                             to indices pindx: \n
01773                             pcoords(i,0) x-coordinate of node i, \n
01774                             pcoords(i,1) y-coordinate of node i.
01775   \param  t         [in]  - Matrix (ELE x 3) of indices of the vertices in a triangle: \n
01776                             t(i,j) index of the j-th vertex in triangle i, where i = 1, ..., ELE
01777   \param  y         [in]  - Reference-counting pointer to the Epetra_MultiVector at which the nonlinear
01778                             term is evaluated.
01779   \param  Gp        [out] - Reference-counting pointer to the Epetra_FECrsMatrix containing the Jacobian
01780                             of the nonlinear form.
01781   \return 0                 if successful.
01782 
01783   \par Detailed Description:
01784 
01785   Assembles the nonlinear term \e Gp, represented by
01786      \f[
01787        \{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,
01788      \f]
01789   where \f$ g'(y_h) \f$ is given in the local function \e gpfctn, and \f$\{ \phi_j \}_{j = 1}^{m}\f$ is the
01790   piecewise linear nodal basis for the state space.
01791 */
01792 int GLpApp::nonlinjac(const Epetra_Comm & Comm,
01793               const Epetra_IntSerialDenseVector & ipindx,
01794               const Epetra_SerialDenseMatrix & ipcoords,
01795               const Epetra_IntSerialDenseVector & pindx,
01796               const Epetra_SerialDenseMatrix & pcoords,
01797               const Epetra_IntSerialDenseMatrix & t,
01798               const Teuchos::RefCountPtr<const Epetra_MultiVector> & y,
01799               Teuchos::RefCountPtr<Epetra_FECrsMatrix> & Gp)
01800 {
01801 
01802   int myPID = Comm.MyPID();
01803   int numProcs = Comm.NumProc();
01804 
01805   int numLocNodes     = pindx.M();
01806   int numMyLocNodes   = ipindx.M();
01807   int numLocElems     = t.M();
01808   int numNodesPerElem = 3;
01809 
01810   int indexBase = 1;
01811 
01812   Epetra_Map standardmap(-1, numMyLocNodes, (int*)ipindx.A(), indexBase, Comm);
01813   Epetra_Map overlapmap(-1, numLocNodes, (int*)pindx.A(), indexBase, Comm);
01814 
01815   int format = Epetra_FECrsMatrix::COLUMN_MAJOR;
01816   Gp = rcp(new Epetra_FECrsMatrix(Copy, standardmap, 0));
01817 
01818   int* nodes = new int[numNodesPerElem];
01819   int i=0, j=0, err=0;
01820   
01821   // get quadrature nodes and weights
01822   Epetra_SerialDenseMatrix Nodes;
01823   Epetra_SerialDenseVector Weights;
01824   quadrature(2,3,Nodes,Weights);
01825   int numQuadPts = Nodes.M();
01826 
01827   // Evaluate nodal basis functions and their derivatives at quadrature points
01828   // N(i,j) = value of the j-th basis function at quadrature node i.
01829   Epetra_SerialDenseMatrix N;
01830   N.Shape(numQuadPts,3);
01831   for (int i=0; i<numQuadPts; i++) {
01832     N(i,0) = 1.0 - Nodes(i,0) - Nodes(i,1);
01833     N(i,1) = Nodes(i,0);
01834     N(i,2) = Nodes(i,1);
01835   }
01836   
01837   // Declare quantities needed for the call to the local assembly routine.
01838   Epetra_IntSerialDenseVector epetra_nodes(View, nodes, numNodesPerElem);
01839   Epetra_SerialDenseMatrix vertices(numNodesPerElem, pcoords.N());
01840 
01841   Epetra_SerialDenseVector ly;          // local entries of y
01842   Epetra_SerialDenseVector Nly;         // N*ly
01843   Epetra_SerialDenseVector lgfctn;      // gfctn(Nly)
01844   Epetra_SerialDenseVector lgfctnNiNj;  // lgfctn.*N(:,i).*N(:,j)
01845   Epetra_SerialDenseMatrix lGp;         // local contribution
01846   // Size and init to zero.
01847   ly.Size(numNodesPerElem);
01848   Nly.Size(numQuadPts);
01849   lgfctn.Size(numQuadPts);
01850   lgfctnNiNj.Size(numQuadPts);
01851   lGp.Shape(numNodesPerElem, numNodesPerElem);
01852   
01853   Epetra_SerialDenseMatrix B(2,2);
01854   double adB;
01855   
01856   for(i=0; i<numLocElems; i++) {
01857 
01858     nodes[0] = t(i,0); nodes[1] = t(i,1); nodes[2] = t(i,2);
01859     for (j=0; j<numNodesPerElem; j++) {
01860       vertices(j,0) = pcoords(overlapmap.LID(nodes[j]), 0);
01861       vertices(j,1) = pcoords(overlapmap.LID(nodes[j]), 1);
01862     }
01863 
01864     // Construct affine transformation matrix.
01865     for(int i=0; i<2; i++) {
01866       B(i,0) = vertices(1,i)-vertices(0,i);
01867       B(i,1) = vertices(2,i)-vertices(0,i);
01868     }
01869     adB  = abs(determinant(B));
01870 
01871     // Construct local (to each processor) element view of y. 
01872     for (j=0; j<numNodesPerElem; j++) {
01873       ly(j) = (*((*y)(0)))[overlapmap.LID(nodes[j])];
01874     }
01875 
01876     Nly.Multiply('N', 'N', 1.0, N, ly, 0.0);
01877     gpfctn(Nly, lgfctn);
01878     
01879     for (int i=0; i<numNodesPerElem; i++) {
01880       for (int j=0; j<numNodesPerElem; j++) {
01881         compproduct(lgfctnNiNj, lgfctn.A(), N[i], N[j]);
01882         lGp(i,j) = adB*lgfctnNiNj.Dot(Weights);
01883       }
01884     }
01885     
01886     err = Gp->InsertGlobalValues(epetra_nodes, lGp, format);
01887     if (err<0) return(err);
01888   }
01889 
01890   // Call global assemble.
01891 
01892   err = Gp->GlobalAssemble();
01893   if (err<0) return(err);
01894 
01895   delete [] nodes;
01896 
01897   return(0);
01898 }
01899 
01900 /*  \brief Componentwise evaluation of the first derivative of the nonlinear reaction term.
01901   \param  v   [in]  - Vector at which the first derivative is evaluated.
01902   \param  gv  [out] - Vector value.
01903 */
01904 void GLpApp::gpfctn(const Epetra_SerialDenseVector & v, Epetra_SerialDenseVector & gv) {
01905   for (int i=0; i<v.M(); i++) {
01906     gv(i) = 3.0*pow(v(i),2)-1.0;
01907   }  
01908 }
01909 
01910 /*  \brief Performs finite-element assembly of the nonlinear reaction term.
01911 
01912   \param  Comm      [in]  - The Epetra (MPI) communicator.
01913   \param  ipindx    [in]  - Vector of NUMIP indices of nodes that are `unique' to a subdomain
01914                             (i.e. owned by the corresponding processor).
01915   \param  ipcoords  [in]  - Matrix (NUMIP x 2) of x-y-coordinates of the vertices cooresponding
01916                             to indices ipindx: \n
01917                             ipcoords(i,0) x-coordinate of node i, \n
01918                             ipcoords(i,1) y-coordinate of node i.
01919   \param  pindx     [in]  - Vector of NUMP indices of ALL nodes in a subdomain, including
01920                             the shared nodes.
01921   \param  pcoords   [in]  - Matrix (NUMP x 2) of x-y-coordinates of the vertices cooresponding
01922                             to indices pindx: \n
01923                             pcoords(i,0) x-coordinate of node i, \n
01924                             pcoords(i,1) y-coordinate of node i.
01925   \param  t         [in]  - Matrix (ELE x 3) of indices of the vertices in a triangle: \n
01926                             t(i,j) index of the j-th vertex in triangle i, where i = 1, ..., ELE
01927   \param  y         [out] - Reference-counting pointer to the Epetra_MultiVector at which the nonlinear
01928                             form is evaluated.
01929   \param  g         [out] - Reference-counting pointer to the Epetra_FEVector containing the value
01930                             of the nonlinear form.
01931   \return 0                 if successful.
01932 
01933   \par Detailed Description:
01934 
01935   Assembles the nonlinear term \e g, represented by
01936      \f[
01937        \{N(y)\}_{j} = \langle g(y_h),\phi_j \rangle =  \int_{\Omega} g(y_h(x)) \phi_j(x) dx,
01938      \f]
01939   where \f$ g(y_h) \f$ is given in the local function \e gfctn, and \f$\{ \phi_j \}_{j = 1}^{m}\f$ is the
01940   piecewise linear nodal basis for the state space.
01941 */
01942 int GLpApp::nonlinvec(const Epetra_Comm & Comm,
01943               const Epetra_IntSerialDenseVector & ipindx,
01944               const Epetra_SerialDenseMatrix & ipcoords,
01945               const Epetra_IntSerialDenseVector & pindx,
01946               const Epetra_SerialDenseMatrix & pcoords,
01947               const Epetra_IntSerialDenseMatrix & t,
01948               const Teuchos::RefCountPtr<const Epetra_MultiVector> & y,
01949               Teuchos::RefCountPtr<Epetra_FEVector> & g)
01950 {
01951 
01952   int myPID = Comm.MyPID();
01953   int numProcs = Comm.NumProc();
01954 
01955   int numLocNodes     = pindx.M();
01956   int numMyLocNodes   = ipindx.M();
01957   int numLocElems     = t.M();
01958   int numNodesPerElem = 3;
01959 
01960   int indexBase = 1;
01961 
01962   Epetra_Map standardmap(-1, numMyLocNodes, (int*)ipindx.A(), indexBase, Comm);
01963   Epetra_Map overlapmap(-1, numLocNodes, (int*)pindx.A(), indexBase, Comm);
01964 
01965   g = rcp(new Epetra_FEVector(standardmap));
01966 
01967   int* nodes = new int[numNodesPerElem];
01968   int i=0, j=0, err=0;
01969   
01970   // get quadrature nodes and weights
01971   Epetra_SerialDenseMatrix Nodes;
01972   Epetra_SerialDenseVector Weights;
01973   quadrature(2,3,Nodes,Weights);
01974   int numQuadPts = Nodes.M();
01975 
01976   // Evaluate nodal basis functions and their derivatives at quadrature points
01977   // N(i,j) = value of the j-th basis function at quadrature node i.
01978   Epetra_SerialDenseMatrix N;
01979   N.Shape(numQuadPts,3);
01980   for (int i=0; i<numQuadPts; i++) {
01981     N(i,0) = 1.0 - Nodes(i,0) - Nodes(i,1);
01982     N(i,1) = Nodes(i,0);
01983     N(i,2) = Nodes(i,1);
01984   }
01985 
01986   // Declare quantities needed for the call to the local assembly routine.
01987   Epetra_IntSerialDenseVector epetra_nodes(View, nodes, numNodesPerElem);
01988   Epetra_SerialDenseMatrix vertices(numNodesPerElem, pcoords.N());
01989 
01990   Epetra_SerialDenseVector ly;        // local entries of y
01991   Epetra_SerialDenseVector Nly;       // N*ly
01992   Epetra_SerialDenseVector lgfctn;    // gfctn(Nly)
01993   Epetra_SerialDenseVector lgfctnNi;  // lgfctn.*N(:,i)
01994   Epetra_SerialDenseVector lg;        // local contribution
01995   // Size and init to zero.
01996   ly.Size(numNodesPerElem);
01997   Nly.Size(numQuadPts);
01998   lgfctn.Size(numQuadPts);
01999   lgfctnNi.Size(numQuadPts);
02000   lg.Size(numNodesPerElem);
02001   
02002   Epetra_SerialDenseMatrix B(2,2);
02003   double adB;
02004   
02005   for(i=0; i<numLocElems; i++) {
02006 
02007     nodes[0] = t(i,0); nodes[1] = t(i,1); nodes[2] = t(i,2);
02008     for (j=0; j<numNodesPerElem; j++) {
02009       vertices(j,0) = pcoords(overlapmap.LID(nodes[j]), 0);
02010       vertices(j,1) = pcoords(overlapmap.LID(nodes[j]), 1);
02011     }
02012 
02013     // Construct affine transformation matrix.
02014     for(int i=0; i<2; i++) {
02015       B(i,0) = vertices(1,i)-vertices(0,i);
02016       B(i,1) = vertices(2,i)-vertices(0,i);
02017     }
02018     adB  = abs(determinant(B));
02019 
02020     // Construct local (to each processor) element view of y. 
02021     for (j=0; j<numNodesPerElem; j++) {
02022       ly(j) = (*((*y)(0)))[overlapmap.LID(nodes[j])];
02023     }
02024 
02025     Nly.Multiply('N', 'N', 1.0, N, ly, 0.0);
02026     gfctn(Nly, lgfctn);
02027     
02028     for (int i=0; i<numNodesPerElem; i++) {
02029       compproduct(lgfctnNi, lgfctn.A(), N[i]);
02030       lg(i) = adB*lgfctnNi.Dot(Weights);
02031     }
02032     
02033     err = g->SumIntoGlobalValues(epetra_nodes, lg);
02034     if (err<0) return(err);
02035   }
02036 
02037   // Call global assemble.
02038 
02039   err = g->GlobalAssemble();
02040   if (err<0) return(err);
02041 
02042   delete [] nodes;
02043 
02044   return(0);
02045 }
02046 
02047 
02048 /*  \brief Componentwise evaluation of the nonlinear reaction term.
02049   \param  v   [in]  - Vector at which the nonlinear function is evaluated.
02050   \param  gv  [out] - Vector value.
02051 */
02052 void GLpApp::gfctn(const Epetra_SerialDenseVector & v, Epetra_SerialDenseVector & gv) {
02053   for (int i=0; i<v.M(); i++) {
02054     gv(i) = pow(v(i),3)-v(i);
02055   }  
02056 }
02057 
02058 /* ======== ================ *
02059  * function CrsMatrix2MATLAB *
02060  * ======== ================ *
02061  *
02062  * Print out a CrsMatrix in a MATLAB format. Each processor prints out
02063  * its part, starting from proc 0 to proc NumProc-1. The first line of
02064  * each processor's output states the number of local rows and of
02065  * local nonzero elements.
02066  *
02067  *
02068  * Return code:        true if matrix has been printed out
02069  * -----------         false otherwise
02070  *
02071  * Parameters:
02072  * ----------
02073  *
02074  * - Epetra_CrsMatrix  reference to the distributed CrsMatrix to 
02075  *                     print out
02076  * - ostream &         reference to output stream 
02077  */
02078 
02079 bool GLpApp::CrsMatrix2MATLAB(const Epetra_CrsMatrix & A, ostream & outfile) 
02080 {
02081 
02082   int MyPID = A.Comm().MyPID(); 
02083   int NumProc = A.Comm().NumProc();
02084 
02085   // work only on transformed matrices;
02086   if( A.IndicesAreLocal() == false ) {
02087     if( MyPID == 0 ) { 
02088       cerr << "ERROR in "<< __FILE__ << ", line " << __LINE__ << endl;
02089       cerr << "Function CrsMatrix2MATLAB accepts\n";
02090       cerr << "transformed matrices ONLY. Please call A.TransformToLoca()\n";
02091       cerr << "on your matrix A to that purpose.\n";
02092       cerr << "Now returning...\n";
02093     }
02094     return false;
02095   }
02096 
02097   int NumMyRows = A.NumMyRows(); // number of rows on this process
02098   int NumNzRow;   // number of nonzero elements for each row
02099   int NumEntries; // number of extracted elements for each row
02100   int NumGlobalRows; // global dimensio of the problem
02101   int GlobalRow;  // row in global ordering
02102   int NumGlobalNonzeros; // global number of nonzero elements
02103 
02104   NumGlobalRows = A.NumGlobalRows();
02105   NumGlobalNonzeros = A.NumGlobalNonzeros();
02106 
02107   // print out on cout if no filename is provided
02108 
02109   int IndexBase = A.IndexBase(); // MATLAB starts from 0
02110   if( IndexBase == 0 )
02111     IndexBase = 1;
02112   else if ( IndexBase == 1)
02113     IndexBase = 0;
02114 
02115   // write on file the dimension of the matrix
02116 
02117   if( MyPID==0 ) {
02118     outfile << "A = spalloc(";
02119     outfile << NumGlobalRows << ',' << NumGlobalRows;
02120     outfile << ',' << NumGlobalNonzeros << ");\n";
02121   }
02122 
02123   for( int Proc=0 ; Proc<NumProc ; ++Proc ) {
02124     A.Comm().Barrier();
02125     if( MyPID == Proc ) {
02126 
02127       outfile << "\n\n% On proc " << Proc << ": ";
02128       outfile << NumMyRows << " rows and ";
02129       outfile << A.NumMyNonzeros() << " nonzeros\n";
02130 
02131       // cycle over all local rows to find out nonzero elements
02132       for( int MyRow=0 ; MyRow<NumMyRows ; ++MyRow ) {
02133 
02134         GlobalRow = A.GRID(MyRow);
02135 
02136         NumNzRow = A.NumMyEntries(MyRow);
02137         double *Values = new double[NumNzRow];
02138         int *Indices = new int[NumNzRow];
02139 
02140         A.ExtractMyRowCopy(MyRow, NumNzRow, 
02141                            NumEntries, Values, Indices);
02142         // print out the elements with MATLAB syntax
02143         for( int j=0 ; j<NumEntries ; ++j ) {
02144           outfile << "A(" << GlobalRow  + IndexBase 
02145                   << "," << A.GCID(Indices[j]) + IndexBase
02146                   << ") = " << Values[j] << ";\n";
02147         }
02148 
02149         delete Values;
02150         delete Indices;
02151       }
02152       
02153     }
02154     A.Comm().Barrier();
02155   }
02156 
02157   return true;
02158 
02159 }
02160 
02161 
02162 /* ======== ============= *
02163  * function Vector2MATLAB *
02164  * ======== ============= *
02165  *
02166  * Print out a Epetra_Vector in a MATLAB format. Each processor prints out
02167  * its part, starting from proc 0 to proc NumProc-1. The first line of
02168  * each processor's output states the number of local rows and of
02169  * local nonzero elements.
02170  *
02171  * Return code:        true if vector has been printed out
02172  * -----------         false otherwise
02173  *
02174  * Parameters:
02175  * ----------
02176  *
02177  * - Epetra_Vector     reference to vector
02178  * - ostream &         reference to output stream 
02179  */
02180 
02181 bool GLpApp::Vector2MATLAB( const Epetra_Vector & v, ostream & outfile)
02182 {
02183   
02184   int MyPID = v.Comm().MyPID(); 
02185   int NumProc = v.Comm().NumProc();
02186   int MyLength = v.MyLength();
02187   int GlobalLength = v.GlobalLength();
02188   
02189   // print out on cout if no filename is provided
02190 
02191   // write on file the dimension of the matrix
02192 
02193   if( MyPID == 0 ) outfile << "v = zeros(" << GlobalLength << ",1)\n";
02194 
02195   int NumMyElements = v.Map().NumMyElements();
02196   // get update list
02197   int * MyGlobalElements = v.Map().MyGlobalElements( );
02198   
02199   int Row;
02200 
02201   int IndexBase = v.Map().IndexBase(); // MATLAB starts from 0
02202   if( IndexBase == 0 )
02203     IndexBase = 1;
02204   else if ( IndexBase == 1)
02205     IndexBase = 0;
02206 
02207   for( int Proc=0 ; Proc<NumProc ; ++Proc ) {
02208     v.Comm().Barrier();
02209     if( MyPID == Proc ) {
02210 
02211       outfile << "% On proc " << Proc << ": ";
02212       outfile << MyLength << " rows of ";
02213       outfile << GlobalLength << " elements\n";
02214 
02215       for( Row=0 ; Row<MyLength ; ++Row ) {
02216         outfile << "v(" << MyGlobalElements[Row] + IndexBase
02217              << ") = " << v[Row] << ";\n";
02218       }
02219       
02220     }
02221       
02222     v.Comm().Barrier();
02223   }
02224 
02225   return true;
02226 
02227 } /* Vector2MATLAB */
02228 
02229 
02230 /* ======== =============== *
02231  * function FEVector2MATLAB *
02232  * ======== =============== *
02233  *
02234  * Print out a Epetra_Vector in a MATLAB format. Each processor prints out
02235  * its part, starting from proc 0 to proc NumProc-1. The first line of
02236  * each processor's output states the number of local rows and of
02237  * local nonzero elements.
02238  *
02239  * Return code:        true if vector has been printed out
02240  * -----------         false otherwise
02241  *
02242  * Parameters:
02243  * ----------
02244  *
02245  * - Epetra_FEVector   reference to FE vector
02246  * - ostream &         reference to output stream 
02247  */
02248 
02249 bool GLpApp::FEVector2MATLAB( const Epetra_FEVector & v, ostream & outfile)
02250 {
02251   
02252   int MyPID = v.Comm().MyPID(); 
02253   int NumProc = v.Comm().NumProc();
02254   int MyLength = v.MyLength();
02255   int GlobalLength = v.GlobalLength();
02256   
02257   // print out on cout if no filename is provided
02258 
02259   // write on file the dimension of the matrix
02260 
02261   if( MyPID == 0 ) outfile << "v = zeros(" << GlobalLength << ",1)\n";
02262 
02263   int NumMyElements = v.Map().NumMyElements();
02264   // get update list
02265   int * MyGlobalElements = v.Map().MyGlobalElements( );
02266   
02267   int Row;
02268 
02269   int IndexBase = v.Map().IndexBase(); // MATLAB starts from 0
02270   if( IndexBase == 0 )
02271     IndexBase = 1;
02272   else if ( IndexBase == 1)
02273     IndexBase = 0;
02274 
02275   for( int Proc=0 ; Proc<NumProc ; ++Proc ) {
02276     v.Comm().Barrier();
02277     if( MyPID == Proc ) {
02278 
02279       outfile << "% On proc " << Proc << ": ";
02280       outfile << MyLength << " rows of ";
02281       outfile << GlobalLength << " elements\n";
02282 
02283       for( Row=0 ; Row<MyLength ; ++Row ) {
02284         outfile << "v(" << MyGlobalElements[Row] + IndexBase
02285              << ") = " << v[0][Row] << ";\n";
02286       }
02287       
02288     }
02289       
02290     v.Comm().Barrier();
02291   }
02292 
02293   return true;
02294 
02295 } /* FEVector2MATLAB */
02296 
02297 
02298 /*  \brief  Returns the nodes and weights for the integration \n
02299              on the interval [0,1] (dim = 1) \n
02300              on the triangle with vertices (0,0), (1,0), (0,1) (if dim = 2) \n
02301              on the tetrahedron with vertices (0,0,0), (1,0,0), (0,1,0), (0,0,1) (if dim = 3).
02302              
02303   \param  dim     [in]   - spatial dimension (dim = 1, 2)
02304   \param  order   [in]   - required degree of polynomials that integrate exactly
02305   \param  nodes   [out]  - Matrix in which the i-th row of nodes gives coordinates of the i-th quadrature node
02306   \param  weights [out]  - quadrature weights
02307 
02308   \return 0                if successful
02309 */
02310 int GLpApp::quadrature(const int dim, const int order,
02311                Epetra_SerialDenseMatrix & nodes,
02312                Epetra_SerialDenseVector & weights)
02313 {
02314   
02315   if (dim == 1) {
02316 
02317     // Gauss quadrature nodes and weights on the interval [0,1]
02318 
02319     if (order == 1) {
02320       nodes.Shape(1,1);
02321       nodes(0,0) = 0.5;
02322       weights.Size(1);
02323       weights(0) = 1.0;
02324     }
02325     else if (order == 2) {
02326       nodes.Shape(2,1);
02327       nodes(0,0) = (1.0-1.0/sqrt(3.0))/2.0;
02328       nodes(1,0) = (1.0+1.0/sqrt(3.0))/2.0;
02329       weights.Size(2);
02330       weights(0) = 0.5;
02331       weights(1) = 0.5;
02332     }
02333     else if (order == 3) {
02334       nodes.Shape(3,1);
02335       nodes(0,0) = (1.0-sqrt(3.0/5.0))/2.0;
02336       nodes(1,0) = 0.5;
02337       nodes(2,0) = (1.0+sqrt(3.0/5.0))/2.0;
02338       weights.Size(3);
02339       weights(0) = 5.0/18.0;
02340       weights(1) = 4.0/9.0;
02341       weights(2) = 5.0/18.0;
02342     }
02343     else {
02344       cout << "Quadrature for dim = " << dim << " and order = ";
02345       cout << order << " not available.\n";
02346       exit(-1);
02347     }
02348 
02349   }
02350   else if (dim == 2) {
02351     
02352     // Quadrature nodes and weights on the unit simplex with
02353     // vertices (0,0), (1,0), and (0,1).
02354 
02355     if (order == 1) {
02356       nodes.Shape(1,2);
02357       nodes(0,0) = 1.0/3.0; nodes (0,1) = 1.0/3.0;
02358       weights.Size(1);
02359       weights(0) = 0.5;
02360     }
02361     else if (order == 2) {
02362       nodes.Shape(3,2);
02363       nodes(0,0) = 1.0/6.0; nodes (0,1) = 1.0/6.0;
02364       nodes(1,0) = 2.0/3.0; nodes (1,1) = 1.0/6.0;
02365       nodes(2,0) = 1.0/6.0; nodes (2,1) = 2.0/3.0;
02366       weights.Size(3);
02367       weights(0) = 1.0/6.0;
02368       weights(1) = 1.0/6.0;
02369       weights(2) = 1.0/6.0;
02370     }
02371     else if (order == 3) {
02372       nodes.Shape(4,2);
02373       nodes(0,0) = 1.0/3.0; nodes (0,1) = 1.0/3.0;
02374       nodes(1,0) = 3.0/5.0; nodes (1,1) = 1.0/5.0;
02375       nodes(2,0) = 1.0/5.0; nodes (2,1) = 3.0/5.0;
02376       nodes(3,0) = 1.0/5.0; nodes (3,1) = 1.0/5.0;
02377       weights.Size(4);
02378       weights(0) = -9.0/32.0;
02379       weights(1) = 25.0/96.0;
02380       weights(2) = 25.0/96.0;
02381       weights(3) = 25.0/96.0;
02382     }
02383     else {
02384       cout << "Quadrature for dim = " << dim << " and order = ";
02385       cout << order << " not available.\n";
02386       exit(-1);
02387     }
02388 
02389   }
02390   else {
02391     cout << "Quadrature for dim = " << dim << " not available.\n";
02392     exit(-1);
02393   }
02394 
02395   return(0);
02396 }

Generated on Thu Sep 18 12:31:44 2008 for EpetraExt by doxygen 1.3.9.1