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

Generated on Wed May 12 21:40:37 2010 for EpetraExt by  doxygen 1.4.7