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

Generated on Wed May 12 21:24:46 2010 for EpetraExt by  doxygen 1.4.7