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