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 ostream& operator <<(ostream &, const Usr_Par &);
00098 
00099 bool CrsMatrix2MATLAB(const Epetra_CrsMatrix &, ostream &);
00100 
00101 bool Vector2MATLAB( const Epetra_Vector &, ostream &);
00102 
00103 bool FEVector2MATLAB( const Epetra_FEVector &, 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       //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       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   bool MakeDataContiguous = false;
00629   EpetraExt::RowMatrix_Transpose transposer( MakeDataContiguous );
00630   Epetra_CrsMatrix & transA = dynamic_cast<Epetra_CrsMatrix&>(transposer(*A_));
00631   Epetra_CrsMatrix & transB = dynamic_cast<Epetra_CrsMatrix&>(transposer(*B_));
00632   Epetra_CrsMatrix & transNpy = dynamic_cast<Epetra_CrsMatrix&>(transposer(*Npy_));
00633   // Insert transpose of A and Npy into Augmat.
00634   for (int i=0; i<nummystates; i++) {
00635     nummyentries = transA.NumMyEntries(i);
00636     values.Resize(nummyentries);
00637     indices.Resize(nummyentries);
00638     transA.ExtractGlobalRowCopy(KKTmapindx.Values()[i], nummyentries, checkentries, values.Values(),
00639                                 indices.Values());
00640     for (int j=0; j<nummyentries; j++)
00641       indices[j] = indices[j]+2*numstates;
00642     if (nummyentries > 0)
00643       Augmat_->InsertGlobalValues(KKTmapindx.Values()[i], nummyentries, values.Values(), 
00644                                   indices.Values());
00645     nummyentries = transNpy.NumMyEntries(i);
00646     values.Resize(nummyentries);
00647     indices.Resize(nummyentries);
00648     transNpy.ExtractGlobalRowCopy(KKTmapindx.Values()[i], nummyentries, checkentries, values.Values(),
00649                                   indices.Values());
00650     for (int j=0; j<nummyentries; j++)
00651       indices[j] = indices[j]+2*numstates;
00652     if (nummyentries > 0)
00653       Augmat_->InsertGlobalValues(KKTmapindx.Values()[i], nummyentries, values.Values(), 
00654                                   indices.Values());
00655   }
00656   // Insert transpose of B into Augmat.
00657   for (int i=0; i<nummystates; i++) {
00658     nummyentries = transB.NumMyEntries(i);
00659     values.Resize(nummyentries);
00660     indices.Resize(nummyentries);
00661     transB.ExtractGlobalRowCopy(KKTmapindx.Values()[i], nummyentries, checkentries, values.Values(),
00662                                 indices.Values());
00663     for (int j=0; j<nummyentries; j++)
00664       indices[j] = indices[j]+2*numstates;
00665     int err = 0;
00666     if (nummyentries > 0)
00667       err = Augmat_->InsertGlobalValues(KKTmapindx.Values()[i]+numstates, nummyentries,
00668                                         values.Values(), indices.Values());
00669     // This will give a nasty message if something goes wrong with the insertion of B transpose.
00670     if (err < 0) {
00671       cout << "Insertion of entries failed:\n";
00672       cout << indices;
00673       cout << nummyentries << endl;
00674       cout << "at row: " << KKTmapindx.Values()[i]+numstates << endl << endl;
00675     }
00676   }
00677 
00678   Augmat_->FillComplete(KKTmap, KKTmap);
00679   // End building the KKT matrix.
00680 
00681 }
00682 
00683 void GLpYUEpetraDataPool::PrintVec( const Teuchos::RefCountPtr<const Epetra_Vector> & x )
00684 {
00685   Vector2MATLAB(*x, cout);
00686 }
00687 
00688 Usr_Par::Usr_Par() {
00689   Epetra_BLAS blas;
00690   Epetra_SerialDenseVector product(4);
00691 
00692   // get nodes and weights
00693   quadrature(2,3,Nodes,Weights);
00694     
00695   // Evaluate nodal basis functions and their derivatives at quadrature
00696   // pts N(i,j) = value of the j-th basis function at quadrature node i.
00697   N.Shape(Nodes.M(),3);
00698   for (int i=0; i<Nodes.M(); i++) {
00699     N(i,0) = 1.0 - Nodes(i,0) - Nodes(i,1);
00700     N(i,1) = Nodes(i,0);
00701     N(i,2) = Nodes(i,1);
00702   }
00703   // Nx1(i,j) partial derrivative of basis function j wrt x1 at node i
00704   Nx1.Shape(Nodes.M(),3);
00705   for (int i=0; i<Nodes.M(); i++) {
00706     Nx1(i,0) = -1.0;
00707     Nx1(i,1) = 1.0;
00708     Nx1(i,2) = 0.0;
00709   }
00710   // Nx2(i,j) partial derrivative of basis function j wrt x2 at node i
00711   Nx2.Shape(Nodes.M(),3);
00712   for (int i=0; i<Nodes.M(); i++) {
00713     Nx2(i,0) = -1.0;
00714     Nx2(i,1) = 0.0;
00715     Nx2(i,2) = 1.0;
00716   }
00717 
00718   // Evaluate integrals of various combinations of partial derivatives
00719   // of the local basis functions (they're constant).
00720   S1.Shape(3,3);
00721   S1(0,0)= 1.0; S1(0,1)=-1.0; S1(0,2)= 0.0;
00722   S1(1,0)=-1.0; S1(1,1)= 1.0; S1(1,2)= 0.0;
00723   S1(2,0)= 0.0; S1(2,1)= 0.0; S1(2,2)= 0.0;
00724   S2.Shape(3,3);
00725   S2(0,0)= 2.0; S2(0,1)=-1.0; S2(0,2)=-1.0;
00726   S2(1,0)=-1.0; S2(1,1)= 0.0; S2(1,2)= 1.0;
00727   S2(2,0)=-1.0; S2(2,1)= 1.0; S2(2,2)= 0.0;
00728   S3.Shape(3,3);
00729   S3(0,0)= 1.0; S3(0,1)= 0.0; S3(0,2)=-1.0;
00730   S3(1,0)= 0.0; S3(1,1)= 0.0; S3(1,2)= 0.0;
00731   S3(2,0)=-1.0; S3(2,1)= 0.0; S3(2,2)= 1.0;
00732     
00733   // Evaluate integrals of basis functions N(i).
00734   Nw.Size(3);
00735   Nw.Multiply('T', 'N', 1.0, N, Weights, 0.0);
00736 
00737   // Evaluate integrals of 9 products N(i)*N(j).
00738   NNw.Shape(3,3);
00739   for (int i=0; i<3; i++) {
00740     for (int j=0; j<3; j++) {
00741       compproduct(product, N[i], N[j]);
00742       NNw(i,j) = blas.DOT(Weights.M(), Weights.A(), product.A());
00743     }
00744   }
00745 
00746   // Evaluate integrals of 27 products N(i)*N(j)*N(k).
00747   NNNw = new Epetra_SerialDenseMatrix[3];
00748   NNNw[0].Shape(3,3); NNNw[1].Shape(3,3); NNNw[2].Shape(3,3); 
00749   for (int i=0; i<3; i++) {
00750     for (int j=0; j<3; j++) {
00751       for (int k=0; k<3; k++) {
00752         compproduct(product, N[i], N[j], N[k]);
00753         NNNw[k](i,j) = blas.DOT(Weights.M(), Weights.A(), product.A());
00754       }
00755     }
00756   }
00757 
00758   // Evaluate integrals of 27 products N(i)*(dN(j)/dx1)*N(k).
00759   NdNdx1Nw = new Epetra_SerialDenseMatrix[3];
00760   NdNdx1Nw[0].Shape(3,3); NdNdx1Nw[1].Shape(3,3); NdNdx1Nw[2].Shape(3,3); 
00761   for (int i=0; i<3; i++) {
00762     for (int j=0; j<3; j++) {
00763       for (int k=0; k<3; k++) {
00764         compproduct(product, N[i], Nx1[j], N[k]);
00765         NdNdx1Nw[k](i,j) = blas.DOT(Weights.M(), Weights.A(), product.A());
00766       }
00767     }
00768   }
00769 
00770   // Evaluate integrals of 27 products N(i)*(dN(j)/dx2)*N(k).
00771   NdNdx2Nw = new Epetra_SerialDenseMatrix[3];
00772   NdNdx2Nw[0].Shape(3,3); NdNdx2Nw[1].Shape(3,3); NdNdx2Nw[2].Shape(3,3); 
00773   for (int i=0; i<3; i++) {
00774     for (int j=0; j<3; j++) {
00775       for (int k=0; k<3; k++) {
00776         compproduct(product, N[i], Nx2[j], N[k]);
00777         NdNdx2Nw[k](i,j) = blas.DOT(Weights.M(), Weights.A(), product.A());
00778       }
00779     }
00780   }
00781 
00782 }
00783 
00784 void Usr_Par::Print(ostream& os) const {
00785   os << endl;
00786   os << "\n\nQuadrature nodes:";
00787   os << "\n-----------------";
00788   Nodes.Print(os);
00789   os << "\n\nQuadrature weights:";
00790   os << "\n-------------------\n";
00791   Weights.Print(os);
00792   os << "\n\nNodal basis functions:";
00793   os << "\n----------------------";
00794   N.Print(os);
00795   os << "\n\nIntegrals of combinations of partial derivatives:";
00796   os << "\n-------------------------------------------------";
00797   S1.Print(os);
00798   S2.Print(os);
00799   S3.Print(os);
00800   os << "\n\nIntegrals of basis functions:";
00801   os << "\n-----------------------------\n";
00802   Nw.Print(os);
00803   os << "\n\nIntegrals of products N(i)*N(j):";
00804   os << "\n--------------------------------\n";
00805   NNw.Print(os);
00806   os << "\n\nIntegrals of products N(i)*N(j)*N(k):";
00807   os << "\n-------------------------------------\n";
00808   NNNw[0].Print(os); NNNw[1].Print(os); NNNw[2].Print(os);
00809   os << "\n\nIntegrals of products N(i)*(dN(j)/dx1)*N(k):";
00810   os << "\n--------------------------------------------\n";
00811   NdNdx1Nw[0].Print(os); NdNdx1Nw[1].Print(os); NdNdx1Nw[2].Print(os);
00812   os << "\n\nIntegrals of products N(i)*(dN(j)/dx2)*N(k):";
00813   os << "\n--------------------------------------------\n";
00814   NdNdx2Nw[0].Print(os); NdNdx2Nw[1].Print(os); NdNdx2Nw[2].Print(os);
00815 }
00816 
00817 ostream& operator <<(ostream & out, const Usr_Par & usr_par)
00818 {
00819   usr_par.Print(out);
00820   return out;
00821 }
00822 
00823 } // namespace GLpApp
00824 
00825 //
00826 // Implementations
00827 //
00828 
00829 int GLpApp::compproduct( Epetra_SerialDenseVector & product,
00830   double *first, double *second)
00831 {
00832   for (int i=0; i<product.M(); i++) {
00833     product[i] = first[i]*second[i];
00834   }
00835   return(0);
00836 }
00837 
00838 int GLpApp::compproduct(Epetra_SerialDenseVector & product,
00839                 double *first, double *second, double *third)
00840 {
00841   for (int i=0; i<product.M(); i++) {
00842     product[i] = first[i]*second[i]*third[i];
00843   }
00844   return(0);
00845 }
00846 
00847 //#define GLPAPP_SHOW_BOUNDARY_ASSEMBLY
00848 
00849 /*  \brief Performs finite-element assembly of face mass matrices.
00850 
00851   \param  Comm      [in]  - The Epetra (MPI) communicator.
00852   \param  ipindx    [in]  - Vector of NUMIP indices of nodes that are `unique' to a subdomain
00853                             (i.e. owned by the corresponding processor).
00854   \param  ipcoords  [in]  - Matrix (NUMIP x 2) of x-y-coordinates of the vertices cooresponding
00855                             to indices ipindx: \n
00856                             ipcoords(i,0) x-coordinate of node i, \n
00857                             ipcoords(i,1) y-coordinate of node i.
00858   \param  pindx     [in]  - Vector of NUMP indices of ALL nodes in a subdomain, including
00859                             the shared nodes.
00860   \param  pcoords   [in]  - Matrix (NUMP x 2) of x-y-coordinates of the vertices cooresponding
00861                             to indices pindx: \n
00862                             pcoords(i,0) x-coordinate of node i, \n
00863                             pcoords(i,1) y-coordinate of node i.
00864   \param  t         [in]  - Matrix (ELE x 3) of indices of the vertices in a triangle: \n
00865                             t(i,j) index of the 0-based j-th vertex in triangle i, where i = 0, ..., numElements-1
00866   \param  e         [in]  - Matrix (EDGE x 3) of edges. \n
00867                             e(i,1:2) contains the indices of the endpoints of edge i, where i = 0, ..., numEdges-1 \n
00868                             e(i,3) contains the boundary marker
00869   \param  B_out     [out] - Reference-counting pointer to the Epetra_FECrsMatrix describing the FE
00870                             state/control face mass matrix.
00871   \param  R_out     [out] - Reference-counting pointer to the Epetra_FECrsMatrix describing the FE
00872                             control/control volume mass matrix.
00873   \return 0                 if successful.
00874 
00875   \par Detailed Description:
00876 
00877   -# Assembles the FE state/control mass matrix \e B, given by
00878      \f[
00879        \mathbf{B}_{jk} = b(\mu_k,\phi_j) = - \int_{\partial \Omega}  \mu_k(x) \phi_j(x) dx,
00880      \f]
00881      where \f$\{ \phi_j \}_{j = 1}^{m}\f$ is the piecewise linear nodal basis for the finite-dimensional
00882      state space, and \f$\{ \mu_j \}_{j = 1}^{n}\f$ is the piecewise linear nodal basis for the
00883      finite-dimensional control space.
00884   -# Assembles the FE control/control mass matrix \e R, given by
00885      \f[
00886        \mathbf{R}_{jk} = b(\mu_k,\mu_j) = - \int_{\partial \Omega}  \mu_k(x) \mu_j(x) dx,
00887      \f]
00888      where \f$\{ \mu_j \}_{j = 1}^{n}\f$ is the piecewise linear nodal basis for the finite-dimensional
00889      control space.
00890 */
00891 int GLpApp::assemble_bdry(
00892   const Epetra_Comm                                &Comm
00893   ,const Epetra_IntSerialDenseVector               &ipindx
00894   ,const Epetra_SerialDenseMatrix                  &ipcoords
00895   ,const Epetra_IntSerialDenseVector               &pindx
00896   ,const Epetra_SerialDenseMatrix                  &pcoords
00897   ,const Epetra_IntSerialDenseMatrix               &t
00898   ,const Epetra_IntSerialDenseMatrix               &e
00899   ,Teuchos::RefCountPtr<Epetra_FECrsMatrix>        *B_out
00900   ,Teuchos::RefCountPtr<Epetra_FECrsMatrix>        *R_out
00901   )
00902 {
00903 
00904   using Teuchos::rcp;
00905 
00906 #ifdef GLPAPP_SHOW_BOUNDARY_ASSEMBLY
00907   Teuchos::RefCountPtr<Teuchos::FancyOStream>
00908     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00909   Teuchos::OSTab tab(out);
00910   *out << "\nEntering assemble_bdry(...) ...\n";
00911 #endif
00912 
00913   int numLocElems = t.M();
00914   int numLocEdges = e.M();
00915 
00916   int indexBase = 1;
00917 
00918   //
00919   // Create a sorted (by global ID) list of boundry nodes
00920   //
00921   int * lastindx = 0;
00922   Epetra_IntSerialDenseVector BdryNode(2*numLocEdges);
00923   for (int j=0; j<numLocEdges; j++) {
00924     BdryNode[j] = e(j,0);
00925     BdryNode[j+numLocEdges] = e(j,1);
00926   }
00927   std::sort(BdryNode.Values(), BdryNode.Values()+2*numLocEdges);
00928   lastindx  = std::unique(BdryNode.Values(), BdryNode.Values()+2*numLocEdges);
00929   const int numBdryNodes = lastindx - BdryNode.Values();
00930   BdryNode.Resize(numBdryNodes); // RAB: This does not overwrite?
00931 
00932   //
00933   // Determine the boundary vertices that belong to this processor.
00934   //
00935   Epetra_IntSerialDenseVector MyBdryNode(numBdryNodes);
00936   lastindx = std::set_intersection(
00937     BdryNode.Values(), BdryNode.Values()+numBdryNodes,  // (Sorted) set A
00938     ipindx.Values(), ipindx.Values()+ipindx.M(),        // (Sorted) set B
00939     MyBdryNode.Values()                                 // Intersection
00940     );
00941   const int numMyBdryNodes = lastindx - MyBdryNode.Values();
00942   MyBdryNode.Resize(numMyBdryNodes); // RAB: This does not overwrite?
00943   
00944   //
00945   // Define the maps for the various lists
00946   //
00947   Epetra_Map standardmap(-1, ipindx.M(), const_cast<int*>(ipindx.A()), indexBase, Comm);
00948   Epetra_Map overlapmap(-1, pindx.M(), const_cast<int*>(pindx.A()), indexBase, Comm);
00949   Epetra_Map mybdryctrlmap(-1, numMyBdryNodes, const_cast<int*>(MyBdryNode.A()), indexBase, Comm);
00950   // Above, it is important to note what mybndyctrlmap represents.  It is the
00951   // a sorted list of global vertex node IDS for nodes on a boundary that are
00952   // uniquely owned by the local process.
00953 
00954 #ifdef GLPAPP_SHOW_BOUNDARY_ASSEMBLY
00955   *out << "\nstandardmap:\n";
00956   standardmap.Print(Teuchos::OSTab(out).o());
00957   *out << "\nmybdryctrlmap:\n";
00958   mybdryctrlmap.Print(Teuchos::OSTab(out).o());
00959 #endif
00960 
00961   //
00962   // Allocate matrices to fill
00963   //
00964   Teuchos::RefCountPtr<Epetra_FECrsMatrix>
00965     B = rcp(new Epetra_FECrsMatrix(Copy,standardmap,0)),
00966     R = rcp(new Epetra_FECrsMatrix(Copy,standardmap,0));
00967   // NOTE: The data map is the same as for the volume matrices. Later, when
00968   // FillComplete is called, we will fix the proper domain and range maps. 
00969   // Declare quantities needed for the call to the local assembly routine.
00970   const int numNodesPerEdge = 2;
00971   Epetra_IntSerialDenseVector epetra_nodes(numNodesPerEdge);
00972 
00973   //
00974   // Load B and R by looping through the edges
00975   //
00976 
00977   Epetra_SerialDenseMatrix Bt(2,2);
00978   int err=0;
00979 
00980   for( int i=0; i < numLocEdges; i++ ) {
00981 
00982     const int
00983       global_id_0 = e(i,0),
00984       global_id_1 = e(i,1),
00985       local_id_0  = overlapmap.LID(global_id_0), // O(log(numip)) bindary search
00986       local_id_1  = overlapmap.LID(global_id_1); // O(log(numip)) bindary search
00987 
00988     epetra_nodes(0) = global_id_0;
00989     epetra_nodes(1) = global_id_1;
00990 
00991     const double
00992       x0 = pcoords(local_id_0,0),
00993       y0 = pcoords(local_id_0,1),
00994       x1 = pcoords(local_id_1,0),
00995       y1 = pcoords(local_id_1,1);
00996     
00997     const double l = sqrt(pow(x0-x1,2) + pow(y0-y1,2));  // Length of this edge
00998     
00999     // We have an explicit formula for Bt.
01000     const double l_sixth = l * (1.0/6.0);
01001     Bt(0,0) = l_sixth * 2.0;
01002     Bt(0,1) = l_sixth * 1.0;
01003     Bt(1,0) = l_sixth * 1.0;
01004     Bt(1,1) = l_sixth * 2.0;
01005 
01006 #ifdef GLPAPP_SHOW_BOUNDARY_ASSEMBLY
01007     *out
01008       << "\nEdge global nodes = ("<<global_id_0<<","<<global_id_1<<"),"
01009       << " local nodes = ("<<local_id_0<<","<<local_id_1<<"),"
01010       << " (x0,y0)(x1,y1)=("<<x0<<","<<y0<<")("<<x1<<","<<y1<<"),"
01011       << " Bt = ["<<Bt(0,0)<<","<<Bt(0,1)<<";"<<Bt(1,0)<<","<<Bt(1,1)<<"]\n";
01012 #endif
01013 
01014     const int format = Epetra_FECrsMatrix::COLUMN_MAJOR;
01015     err = B->InsertGlobalValues(epetra_nodes,Bt,format);
01016     if (err<0) return(err);
01017     err = R->InsertGlobalValues(epetra_nodes,Bt,format);
01018     if (err<0) return(err);
01019     
01020   }
01021 
01022 /*
01023 
01024   err = B->GlobalAssemble(false);
01025   if (err<0) return(err);
01026   err = R->GlobalAssemble(false);
01027   if (err<0) return(err);
01028 
01029   err = B->FillComplete(mybdryctrlmap,standardmap);
01030   if (err<0) return(err);
01031   err = R->FillComplete(mybdryctrlmap,mybdryctrlmap);
01032   if (err<0) return(err);
01033 
01034 */
01035 
01036   err = B->GlobalAssemble(mybdryctrlmap,standardmap,true);
01037   if (err<0) return(err);
01038   err = R->GlobalAssemble(mybdryctrlmap,mybdryctrlmap,true);
01039   if (err<0) return(err);
01040 
01041   if(B_out) *B_out = B;
01042   if(R_out) *R_out = R;
01043 
01044 #ifdef GLPAPP_SHOW_BOUNDARY_ASSEMBLY
01045   *out << "\nB =\n";
01046   B->Print(Teuchos::OSTab(out).o());
01047   *out << "\nLeaving assemble_bdry(...) ...\n";
01048 #endif
01049 
01050   return(0);
01051 
01052 }
01053 
01054 /*  \brief Performs finite-element assembly of volume stiffness and mass matrices,
01055             and the right-hand side (forcing term).
01056 
01057   \param  Comm      [in]  - The Epetra (MPI) communicator.
01058   \param  ipindx    [in]  - Vector of NUMIP indices of nodes that are `unique' to a subdomain
01059                             (i.e. owned by the corresponding processor).
01060   \param  ipcoords  [in]  - Matrix (NUMIP x 2) of x-y-coordinates of the vertices cooresponding
01061                             to indices ipindx: \n
01062                             ipcoords(i,0) x-coordinate of node i, \n
01063                             ipcoords(i,1) y-coordinate of node i.
01064   \param  pindx     [in]  - Vector of NUMP indices of ALL nodes in a subdomain, including
01065                             the shared nodes.
01066   \param  pcoords   [in]  - Matrix (NUMP x 2) of x-y-coordinates of the vertices cooresponding
01067                             to indices pindx: \n
01068                             pcoords(i,0) x-coordinate of node i, \n
01069                             pcoords(i,1) y-coordinate of node i.
01070   \param  t         [in]  - Matrix (ELE x 3) of indices of the vertices in a triangle: \n
01071                             t(i,j) index of the j-th vertex in triangle i, where i = 1, ..., ELE
01072   \param  e         [in]  - Matrix (EDGE x 3) of edges. \n
01073                             e(i,1:2) contains the indices of the endpoints of edge i, where i = 1, ..., EDGE \n
01074                             e(i,3) contains the boundary marker \n
01075                             e(i,3) = 1  Dirichlet boundary conditions on edge i \n
01076                             e(i,3) = 2  Neumann boundary conditions on edge i \n
01077                             e(i,3) = 3  Robin boundary conditions on edge i
01078   \param  A         [out] - Reference-counting pointer to the Epetra_FECrsMatrix describing the FE volume
01079                             stiffness matrix for the advection-diffusion equation. Includes advection,
01080                             diffusion, and reaction terms, and modifications that account for the boundary
01081                             conditions.
01082   \param  M         [out] - Reference-counting pointer to the Epetra_FECrsMatrix describing the FE volume
01083                             mass matrix.
01084   \param  b         [out] - Reference-counting pointer to the Epetra_FEVector describing the discretized
01085                             forcing term in the advection-diffusion equation. Includes modifications that
01086                             account for the boundary conditions.
01087   \return 0                 if successful.
01088 
01089   \par Detailed Description:
01090 
01091   -# Assembles the FE volume stiffness matrix \e A and the right-hand side \e b for the
01092   solution of an advection-diffusion equation using piecewise linear finite elements.
01093   The advection-diffusion equation is 
01094   \f{align*}
01095        - \nabla \cdot \left( K(x) \nabla y(x) \right) + \mathbf{c}(x) \cdot \nabla y(x) + r(x) y(x)
01096        &= f(x), & x &\in \Omega,\;  \\
01097        y(x) &= d(x),    & x &\in {\partial \Omega}_d, \\
01098        K(x) \frac{\partial}{\partial \mathbf{n}}  y(x) &= g(x), & x &\in {\partial \Omega}_n, \\
01099        \sigma_0 K(x) \frac{\partial}{\partial \mathbf{n}} y(x)
01100        + \sigma_1 y(x) &= g(x), & x &\in {\partial \Omega}_r,
01101   \f}
01102   where \f$ K \f$ represents scalar diffusion, \f$ \mathbf{c} \f$ is the advection vector field,
01103   \f$ r \f$ is the reaction term, \f$ d \f$ and  \f$ g \f$ are given functions, \f$sigma_0\f$ and
01104   \f$ sigma_1 \f$ are given quantities, \f$ {\partial \Omega}_d \f$ is the Dirichlet boundary,
01105   \f$ {\partial \Omega}_n \f$ is the Neumann boundary, and \f$ {\partial \Omega}_r \f$ is the Robin
01106   boundary. The quantities \f$ K \f$, \f$ \mathbf{c} \f$, \f$ r \f$, \f$ d \f$, and \f$ g \f$ are
01107   assumed to be piecewise linear. Currently, they are to be hard-coded inside this function.
01108   -# Assembles the FE volume mass matrix \e M.
01109 */
01110 int GLpApp::assemble(const Epetra_Comm & Comm,
01111              const Epetra_IntSerialDenseVector & ipindx,
01112              const Epetra_SerialDenseMatrix & ipcoords,
01113              const Epetra_IntSerialDenseVector & pindx,
01114              const Epetra_SerialDenseMatrix & pcoords,
01115              const Epetra_IntSerialDenseMatrix & t,
01116              const Epetra_IntSerialDenseMatrix & e,
01117              RefCountPtr<Epetra_FECrsMatrix> & A,
01118              RefCountPtr<Epetra_FECrsMatrix> & M,
01119              RefCountPtr<Epetra_FEVector> & b)
01120 {
01121 
01122   int myPID = Comm.MyPID();
01123   int numProcs = Comm.NumProc();
01124   Usr_Par usr_par;
01125 
01126   int numLocElems = t.M();
01127   int numNodesPerElem = 3;
01128 
01129   int indexBase = 1;
01130 
01131   Epetra_Map standardmap(-1, ipindx.M(), (int*)ipindx.A(), indexBase, Comm);
01132   Epetra_Map overlapmap(-1, pindx.M(), (int*)pindx.A(), indexBase, Comm);
01133 
01134   int* nodes = new int[numNodesPerElem];
01135   int i=0, j=0, err=0;
01136 
01137   A = rcp(new Epetra_FECrsMatrix(Copy, standardmap, 0));
01138   M = rcp(new Epetra_FECrsMatrix(Copy, standardmap, 0));
01139   b = rcp(new Epetra_FEVector(standardmap,1));
01140 
01141   // Declare quantities needed for the call to the local assembly routine.
01142   int format = Epetra_FECrsMatrix::COLUMN_MAJOR;
01143   Epetra_IntSerialDenseVector epetra_nodes(View, nodes, numNodesPerElem);
01144 
01145 
01146   /* ************************  First, we build A and b.  ************************ */
01147   Epetra_SerialDenseMatrix At;
01148   Epetra_SerialDenseVector bt;
01149   Epetra_SerialDenseMatrix vertices(numNodesPerElem, pcoords.N());
01150   
01151   Epetra_SerialDenseVector k(numNodesPerElem);
01152   for (i=0; i< numNodesPerElem; i++) k(i)=1.0;
01153   Epetra_SerialDenseMatrix c(numNodesPerElem,2);
01154   for (i=0; i< numNodesPerElem; i++) { c(i,0)=0.0; c(i,1)=0.0; }
01155   Epetra_SerialDenseVector r(numNodesPerElem);
01156   for (i=0; i< numNodesPerElem; i++) r(i)=0.0;
01157   Epetra_SerialDenseVector f(numNodesPerElem);
01158   for (i=0; i< numNodesPerElem; i++) f(i)=0.0;
01159   Epetra_SerialDenseVector g(2);
01160   g(0)=0.0; g(1)=0.0;
01161   Epetra_SerialDenseVector sigma(2);
01162   sigma(0)=0.0; sigma(1)=0.0;
01163   for(i=0; i<numLocElems; i++) {
01164     nodes[0] = t(i,0); nodes[1] = t(i,1); nodes[2] = t(i,2);
01165     for (j=0; j<numNodesPerElem; j++) {
01166       // get vertex
01167       vertices(j,0) = pcoords(overlapmap.LID(nodes[j]), 0);
01168       vertices(j,1) = pcoords(overlapmap.LID(nodes[j]), 1);
01169       // set rhs function
01170       f(j) = cos(GLp_pi*vertices(j,0))*cos(GLp_pi*vertices(j,1)) * 
01171                (2*GLp_pi*GLp_pi + pow(cos(GLp_pi*vertices(j,0)),2)*pow(cos(GLp_pi*vertices(j,1)),2) - 1.0);
01172     }
01173     lassembly(vertices, k, c, r, f, usr_par, At, bt);
01174     err = A->InsertGlobalValues(epetra_nodes, At, format);
01175     if (err<0) return(err);
01176     err = b->SumIntoGlobalValues(epetra_nodes, bt);
01177     if (err<0) return(err);
01178   }
01179 
01180   // MAKE ADJUSTMENTS TO A and b TO ACCOUNT FOR BOUNDARY CONDITIONS.
01181 
01182   // Get Neumann boundary edges.
01183   Epetra_IntSerialDenseMatrix neumann(e.M(),2);
01184   j = 0;
01185   for (i=0; i<e.M(); i++) {
01186     if (e(i,2)==2) {
01187       neumann(j,0) = e(i,0);  neumann(j,1) = e(i,1);
01188       j++;
01189     }
01190   }
01191   neumann.Reshape(j,2);
01192   // Adjust for Neumann BC's.
01193   if (neumann.M() != 0) {
01194     Epetra_BLAS blas;
01195     Epetra_SerialDenseMatrix quadnodes;
01196     Epetra_SerialDenseVector quadweights;
01197     Epetra_SerialDenseMatrix N;
01198     Epetra_SerialDenseMatrix NN;
01199     Epetra_SerialDenseVector product(2);
01200 
01201     // Get quadrature weights and points.
01202     quadrature(1, 2, quadnodes, quadweights);
01203     
01204     // Evaluate nodal basis functions at quadrature points
01205     // N(i,j) value of basis function j at quadrature node i
01206     N.Shape(quadnodes.M(),2);
01207     for (i=0; i<quadnodes.M(); i++) {
01208       N(i,0) = 1.0 - quadnodes(i,0);
01209       N(i,1) = quadnodes(i,0);
01210     }
01211 
01212     // Evaluate integrals of 4 products N(i)*N(j).
01213     NN.Shape(2,2);
01214     for (i=0; i<2; i++) {
01215       for (j=0; j<2; j++) {
01216         compproduct(product, N[i], N[j]);
01217         NN(i,j) = blas.DOT(quadweights.M(), quadweights.A(), product.A());
01218       }
01219     }
01220 
01221     Epetra_IntSerialDenseVector neumann_nodes(2);
01222     Epetra_SerialDenseVector neumann_add(2);
01223     double l;
01224     for (i=0; i<neumann.M(); i++) {
01225       neumann_nodes(0) = neumann(i,0); neumann_nodes(1) = neumann(i,1);
01226       neumann_add(0) = pcoords(overlapmap.LID(neumann_nodes(0)),0)
01227                       -pcoords(overlapmap.LID(neumann_nodes(1)),0);
01228       neumann_add(1) = pcoords(overlapmap.LID(neumann_nodes(0)),1)
01229                       -pcoords(overlapmap.LID(neumann_nodes(1)),1);
01230       l = blas.NRM2(neumann_add.M(), neumann_add.A());
01231       neumann_add.Multiply('N', 'N', 1.0, NN, g, 0.0);
01232       neumann_add.Scale(l);
01233       err = b->SumIntoGlobalValues(neumann_nodes, neumann_add); 
01234       if (err<0) return(err);
01235     }
01236   }
01237 
01238   // Get Robin boundary edges.
01239   Epetra_IntSerialDenseMatrix robin(e.M(),2);
01240   j = 0;
01241   for (i=0; i<e.M(); i++) {
01242     if (e(i,2)==3) {
01243       robin(j,0) = e(i,0);  robin(j,1) = e(i,1);
01244       j++;
01245     }
01246   }
01247   robin.Reshape(j,2);
01248   // Adjust for Robin BC's.
01249   if (robin.M() != 0) {
01250     Epetra_BLAS blas;
01251     Epetra_SerialDenseMatrix quadnodes;
01252     Epetra_SerialDenseVector quadweights;
01253     Epetra_SerialDenseMatrix N;
01254     Epetra_SerialDenseMatrix NN;
01255     Epetra_SerialDenseMatrix * NNN;
01256     Epetra_SerialDenseVector product(2);
01257 
01258     // Get quadrature weights and points.
01259     quadrature(1, 2, quadnodes, quadweights);
01260     
01261     // Evaluate nodal basis functions at quadrature points
01262     // N(i,j) value of basis function j at quadrature node i
01263     N.Shape(quadnodes.M(),2);
01264     for (i=0; i<quadnodes.M(); i++) {
01265       N(i,0) = 1.0 - quadnodes(i,0);
01266       N(i,1) = quadnodes(i,0);
01267     }
01268 
01269     // Evaluate integrals of 4 products N(i)*N(j).
01270     NN.Shape(2,2);
01271     for (i=0; i<2; i++) {
01272       for (j=0; j<2; j++) {
01273         compproduct(product, N[i], N[j]);
01274         NN(i,j) = blas.DOT(quadweights.M(), quadweights.A(), product.A());
01275       }
01276     }
01277 
01278     // Evaluate integrals of 8 products N(i)*N(j)*N(k).
01279     NNN = new Epetra_SerialDenseMatrix[2];
01280     NNN[0].Shape(2,2); NNN[1].Shape(2,2);
01281     for (i=0; i<2; i++) {
01282       for (j=0; j<2; j++) {
01283         for (int k=0; k<2; k++) {
01284           compproduct(product, N[i], N[j], N[k]);
01285           NNN[k](i,j) = blas.DOT(quadweights.M(), quadweights.A(),
01286                                  product.A());
01287         }
01288       }
01289     }
01290 
01291     Epetra_IntSerialDenseVector robin_nodes(2);
01292     Epetra_SerialDenseVector robin_b_add(2);
01293     Epetra_SerialDenseMatrix robin_A_add(2,2);
01294     double l;
01295     for (i=0; i<robin.M(); i++) {
01296       robin_nodes(0) = robin(i,0); robin_nodes(1) = robin(i,1);
01297       
01298       robin_b_add(0) = pcoords(overlapmap.LID(robin_nodes(0)),0)
01299                       -pcoords(overlapmap.LID(robin_nodes(1)),0);
01300       robin_b_add(1) = pcoords(overlapmap.LID(robin_nodes(0)),1)
01301                       -pcoords(overlapmap.LID(robin_nodes(1)),1);
01302       l = blas.NRM2(robin_b_add.M(), robin_b_add.A());
01303       robin_b_add.Multiply('N', 'N', 1.0, NN, g, 0.0);
01304       robin_b_add.Scale(l);
01305       err = b->SumIntoGlobalValues(robin_nodes, robin_b_add); 
01306       if (err<0) return(err);
01307 
01308       NNN[0].Scale(sigma(0)); NNN[1].Scale(sigma(1));
01309       robin_A_add += NNN[0]; robin_A_add += NNN[1];
01310       robin_A_add.Scale(l);
01311       err = A->InsertGlobalValues(robin_nodes, robin_A_add, format);
01312       if (err<0) return(err);
01313     }
01314 
01315     delete [] NNN;
01316   }
01317 
01318   // Get Dirichlet boundary edges.
01319   Epetra_IntSerialDenseMatrix dirichlet(e.M(),2);
01320   j = 0;
01321   for (i=0; i<e.M(); i++) {
01322     if (e(i,2)==1) {
01323       dirichlet(j,0) = e(i,0);  dirichlet(j,1) = e(i,1);
01324       j++;
01325     }
01326   }
01327   dirichlet.Reshape(j,2);
01328   // DIRICHLET NOT DONE! DO THIS LATER!!!!
01329 
01330   /* ************************  Done building A and b.  ************************ */
01331 
01332 
01333 
01334   /* ************************  Second, we build M.  ************************ */
01335 
01336   Epetra_SerialDenseMatrix Mt;
01337 
01338   for (i=0; i< numNodesPerElem; i++) k(i)=0.0;
01339   for (i=0; i< numNodesPerElem; i++) { c(i,0)=0.0; c(i,1)=0.0; }
01340   for (i=0; i< numNodesPerElem; i++) r(i)=1.0;
01341   for (i=0; i< numNodesPerElem; i++) f(i)=0.0;
01342   g(0)=0.0; g(1)=0.0;
01343   sigma(0)=0.0; sigma(1)=0.0;
01344   for(i=0; i<numLocElems; i++) {
01345     nodes[0] = t(i,0); nodes[1] = t(i,1); nodes[2] = t(i,2);
01346     for (j=0; j<numNodesPerElem; j++) {
01347       vertices(j,0) = pcoords(overlapmap.LID(nodes[j]), 0);
01348       vertices(j,1) = pcoords(overlapmap.LID(nodes[j]), 1);
01349     }
01350     lassembly(vertices, k, c, r, f, usr_par, Mt, bt);
01351     err = M->InsertGlobalValues(epetra_nodes, Mt, format);
01352     if (err<0) return(err);
01353   }
01354 
01355   /* ************************  Done building M.  ************************ */
01356 
01357 
01358 
01359   // Call global assemble and FillComplete at the same time.
01360 
01361   err = A->GlobalAssemble();
01362   if (err<0) return(err);
01363 
01364   err = b->GlobalAssemble();
01365   if (err<0) return(err);
01366 
01367   err = M->GlobalAssemble();
01368   if (err<0) return(err);
01369 
01370   delete [] nodes;
01371 
01372   return(0);
01373 }
01374 
01375 /* \brief Computes local stiffness matrix and local RHS vector for simplex
01376            (triangles in two dimensions).
01377                             
01378   \param  vertices   [in]  - Matrix containing the global coordinates of the current simplex.
01379   \param  k          [in]  - Vector containing the value of the diffusion \f$k\f$ at each vertex
01380                              (\f$k\f$ is assumed to be piecewise linear), where \f$k\f$ is the coeff in the
01381                              term \f$ \nabla \cdot (k \nabla(u)) \f$ in the advection-diffusion equation.
01382   \param  c          [in]  - Matrix containing the value of the advection field \f$ \mathbf{c} \f$ at each
01383                              vertex (\f$ \mathbf{c} \f$ is assumed to be piecewise linear), where
01384                              \f$ \mathbf{c} \f$ is the 2-vector of coeffs in the term
01385                              \f$ \mathbf{c}(x) \cdot \nabla y(x) \f$ in the advection-diffusion equation.
01386   \param  r          [in]  - Vector containing the value of \f$ r \f$ at each vertex (\f$ r \f$ is assumed
01387                              to be piecewise linear), where \f$ r \f$ is the coefficient in the term
01388                              \f$ ru \f$ in the advection-diffusion equation.
01389   \param  f          [in]  - Vector containing the value of \f$ f \f$ at each vertex (\f$ f \f$ is assumed to be
01390                              piecewise linear), where \f$ f \f$ is the right hand side in the adv-diff eq
01391   \param  usr_par    [in]  - class containing: \n
01392                               - S1, S2, S3 (3x3) the integrals of various combinations of partials
01393                                 of local basis functions
01394                               - N (1x3) integrals of local basis functions
01395                               - NNN[3] (3x3), integrals of products of local basis functions N(i)*N(j)*N(k)
01396                               - etc.
01397   \param  At         [out] - stiffness matrix contribution for the simplex
01398   \param  bt         [out] - right-hand side contribution for the simplex
01399 
01400   \return 0                 if successful.
01401 
01402   \par Detailed Description:
01403 
01404      Computes the local (per simplex) contributions to the FE volume stiffness matrix \e A and the
01405      right-hand side \e b for the advection-diffusion equation
01406     \f{align*}
01407        - \nabla \cdot \left( K(x) \nabla y(x) \right) + \mathbf{c}(x) \cdot \nabla y(x) + r(x) y(x)
01408        &= f(x), & x &\in \Omega,\;  \\
01409        y(x) &= d(x),    & x &\in {\partial \Omega}_d, \\
01410        K(x) \frac{\partial}{\partial \mathbf{n}}  y(x) &= g(x), & x &\in {\partial \Omega}_n, \\
01411        \sigma_0 K(x) \frac{\partial}{\partial \mathbf{n}} y(x)
01412        + \sigma_1 y(x) &= g(x), & x &\in {\partial \Omega}_r.
01413     \f}
01414 */
01415 int GLpApp::lassembly(const Epetra_SerialDenseMatrix & vertices,
01416               const Epetra_SerialDenseVector & k,
01417               const Epetra_SerialDenseMatrix & c,
01418               const Epetra_SerialDenseVector & r,
01419               const Epetra_SerialDenseVector & f,
01420               const Usr_Par & usr_par,
01421               Epetra_SerialDenseMatrix & At,
01422               Epetra_SerialDenseVector & bt)
01423 {
01424   // Note that the constructors below initialize all entries to 0.
01425   Epetra_SerialDenseMatrix B(2,2);
01426   Epetra_SerialDenseMatrix b(2,2);
01427   Epetra_SerialDenseMatrix BtB(2,2);  
01428   Epetra_SerialDenseMatrix C(2,2);  
01429   Epetra_SerialDenseMatrix M1(3,3);
01430   Epetra_SerialDenseMatrix M2(3,3);
01431   Epetra_SerialDenseMatrix M3(3,3);
01432   Epetra_SerialDenseMatrix tmp(3,3);
01433   double dB, adB;
01434   At.Shape(3,3);
01435   bt.Size(3);
01436 
01437   // Construct affine transformation matrix.
01438   for(int i=0; i<2; i++) {
01439     B(i,0) = vertices(1,i)-vertices(0,i);
01440     B(i,1) = vertices(2,i)-vertices(0,i);
01441   }
01442   dB  = determinant(B);
01443   adB = abs(dB);
01444 
01445   // Compute matrix C = inv(B'*B).
01446   BtB.Multiply('T', 'N', 1.0, B, B, 0.0);
01447   inverse(BtB, C);
01448 
01449   inverse(B, b); b.Scale(dB);
01450   
01451   // Compute the part corresponding to div(K*grad(u)).
01452   tmp = usr_par.S1; tmp.Scale(C(0,0));
01453   M1 += tmp;
01454   tmp = usr_par.S2; tmp.Scale(C(0,1));
01455   M1 += tmp;
01456   tmp = usr_par.S3; tmp.Scale(C(1,1));
01457   M1 += tmp;
01458   M1.Scale( (k(0)*usr_par.Nw(0) + k(1)*usr_par.Nw(1) +
01459              k(2)*usr_par.Nw(2)) * adB );
01460 
01461   // Compute the part corresponding to c'*grad(u).
01462   tmp = usr_par.NdNdx1Nw[0]; tmp.Scale(b(0,0)*c(0,0)+b(0,1)*c(0,1));
01463   M2 += tmp;
01464   tmp = usr_par.NdNdx1Nw[1]; tmp.Scale(b(0,0)*c(1,0)+b(0,1)*c(1,1));
01465   M2 += tmp;
01466   tmp = usr_par.NdNdx1Nw[2]; tmp.Scale(b(0,0)*c(2,0)+b(0,1)*c(2,1));
01467   M2 += tmp;
01468   tmp = usr_par.NdNdx2Nw[0]; tmp.Scale(b(1,0)*c(0,0)+b(1,1)*c(0,1));
01469   M2 += tmp;
01470   tmp = usr_par.NdNdx2Nw[1]; tmp.Scale(b(1,0)*c(1,0)+b(1,1)*c(1,1));
01471   M2 += tmp;
01472   tmp = usr_par.NdNdx2Nw[2]; tmp.Scale(b(1,0)*c(2,0)+b(1,1)*c(2,1));
01473   M2 += tmp;
01474   M2.Scale(adB/dB);
01475 
01476   // Compute the part corresponding to r*u.
01477   tmp = usr_par.NNNw[0]; tmp.Scale(r(0));
01478   M3 += tmp;
01479   tmp = usr_par.NNNw[1]; tmp.Scale(r(1));
01480   M3 += tmp;
01481   tmp = usr_par.NNNw[2]; tmp.Scale(r(2));
01482   M3 += tmp;
01483   M3.Scale(adB);
01484 
01485   At += M1;
01486   At += M2;
01487   At += M3;
01488 
01489   bt.Multiply('T', 'N', adB, usr_par.NNw, f, 0.0);  
01490   
01491   return(0);
01492 }
01493 
01494 /*  \brief Computes the inverse of a dense matrix.
01495 
01496   \param  mat  [in]  - the matrix that is to be inverted
01497   \param  inv  [in]  - the inverse
01498 
01499   \return 0            if successful
01500 */
01501 int GLpApp::inverse(const Epetra_SerialDenseMatrix & mat,
01502             Epetra_SerialDenseMatrix & inv)
01503 {
01504   Epetra_LAPACK lapack;
01505   int dim = mat.M();
01506   int info;
01507   Epetra_IntSerialDenseVector ipiv(dim);
01508   Epetra_SerialDenseVector work(2*dim);
01509 
01510   inv.Shape(dim,dim);
01511   inv = mat;
01512 
01513   lapack.GETRF(dim, dim, inv.A(), dim, ipiv.A(), &info);
01514   lapack.GETRI(dim, inv.A(), dim, ipiv.A(), work.A(), &dim, &info);
01515   
01516   return(0);
01517 }
01518 
01519 /*  \brief Computes the determinant of a dense matrix.
01520 
01521   \param  mat  [in]  - the matrix
01522 
01523   \return the determinant
01524 */
01525 double GLpApp::determinant(const Epetra_SerialDenseMatrix & mat)
01526 {
01527   //Teuchos::LAPACK<int,double> lapack;
01528   Epetra_LAPACK lapack;
01529   double det;
01530   int swaps; 
01531   int dim = mat.M();
01532   int info;
01533   Epetra_IntSerialDenseVector ipiv(dim);
01534  
01535   Epetra_SerialDenseMatrix mymat(mat);
01536   lapack.GETRF(dim, dim, mymat.A(), dim, ipiv.A(), &info);
01537 
01538   det = 1.0;
01539   for (int i=0; i<dim; i++) {
01540     det *= mymat(i,i);
01541   }
01542 
01543   swaps = 0;
01544   for (int i=0; i<dim; i++) {
01545     if ((ipiv[i]-1) != i)
01546       swaps++;
01547   }
01548 
01549   det *= pow((double)(-1.0),swaps);
01550 
01551   return(det);
01552 }
01553 
01554 int GLpApp::meshreader(const Epetra_Comm & Comm,
01555                Epetra_IntSerialDenseVector & ipindx,
01556                Epetra_SerialDenseMatrix & ipcoords,
01557                Epetra_IntSerialDenseVector & pindx,
01558                Epetra_SerialDenseMatrix & pcoords,
01559                Epetra_IntSerialDenseMatrix & t,
01560                Epetra_IntSerialDenseMatrix & e,
01561                const char geomFileBase[],
01562                const bool trace,
01563                const bool dumpAll
01564                )
01565 {
01566   int MyPID = Comm.MyPID();
01567 
01568   const int FileNameSize = 120;
01569   char FileName[FileNameSize];
01570   TEUCHOS_TEST_FOR_EXCEPT(static_cast<int>(std::strlen(geomFileBase) + 5) > FileNameSize);
01571   sprintf(FileName, "%s.%03d", geomFileBase, MyPID);
01572 
01573   {
01574     std::ifstream file_in(FileName);
01575     TEUCHOS_TEST_FOR_EXCEPTION(
01576       file_in.eof(), std::logic_error
01577       ,"Error, the file \""<<FileName<<"\" could not be opened for input!"
01578       );
01579   }
01580 
01581   FILE* fid;
01582   fid = fopen(FileName, "r");
01583 
01584   if(trace) printf("\nReading node info from %s ...\n", FileName);
01585   int numip = 0, numcp = 0; // # owned nodes, # shared nodes
01586   fscanf(fid, "%d %d", &numip, &numcp);
01587   const int nump = numip + numcp; // # total nodes
01588   if(trace) printf( "\nnumip = %d, numcp = %d, nump = %d\n", numip, numcp, nump );
01589   ipindx.Size(numip);
01590   ipcoords.Shape(numip, 2);
01591   pindx.Size(nump);
01592   pcoords.Shape(nump, 2);
01593   for (int i=0; i<numip; i++) {
01594     fscanf(fid,"%d %lf %lf %*d",&ipindx(i),&ipcoords(i,0),&ipcoords(i,1));
01595     if(trace&&dumpAll) printf("%d %lf %lf\n",ipindx(i),ipcoords(i,0),ipcoords(i,1));
01596     pindx(i) = ipindx(i);
01597     pcoords(i,0) = ipcoords(i,0); pcoords(i,1) = ipcoords(i,1);
01598   }
01599   for (int i=numip; i<nump; i++) {
01600     fscanf(fid,"%d %lf %lf %*d",&pindx(i),&pcoords(i,0),&pcoords(i,1));
01601   }
01602 
01603   fscanf(fid,"%*[^\n]");   // Skip to the End of the Line
01604   fscanf(fid,"%*1[\n]");   // Skip One Newline
01605 
01606   fscanf(fid,"%*[^\n]");   // Skip to the End of the Line
01607   fscanf(fid,"%*1[\n]");   // Skip One Newline
01608 
01609   for (int i=0; i<nump; i++) {
01610     fscanf(fid,"%*[^\n]"); // Skip to the End of the Line
01611     fscanf(fid,"%*1[\n]"); // Skip One Newline
01612   }
01613 
01614   if(trace) printf("\nReading element info from %s ...\n", FileName);
01615   int numelems = 0; // # elements on this processor
01616   fscanf(fid, "%d", &numelems);
01617   if(trace) printf( "\nnumelems = %d\n", numelems );
01618   t.Shape(numelems,3);
01619   for (int i=0; i<numelems; i++) {
01620     fscanf(fid, "%d %d %d", &t(i,0), &t(i,1), &t(i,2));
01621     if(trace&&dumpAll) printf("%d %d %d\n", t(i,0), t(i,1), t(i,2));
01622   }
01623 
01624   if(trace) printf("\nReading boundry edge info from %s ...\n", FileName);
01625   int numedges = 0; // # boundry edges on this processor
01626   fscanf(fid,"%d",&numedges);
01627   if(trace) printf( "\nnumedges = %d\n", numedges );
01628   e.Shape(numedges,3);
01629   for (int i=0; i<numedges; i++) {
01630     fscanf(fid, "%d %d %d", &e(i,0), &e(i,1), &e(i,2));
01631     if(trace&&dumpAll) printf("%d %d %d\n", e(i,0), e(i,1), e(i,2));
01632   }
01633 
01634   fclose(fid);
01635   if(trace) printf("Done reading mesh.\n");
01636 
01637   return(0);
01638 
01639 }
01640 
01641 /*  \brief Performs finite-element assembly of the Hessian of the nonlinear form times an arbitrary vector.
01642 
01643   \param  Comm      [in]  - The Epetra (MPI) communicator.
01644   \param  ipindx    [in]  - Vector of NUMIP indices of nodes that are `unique' to a subdomain
01645                             (i.e. owned by the corresponding processor).
01646   \param  ipcoords  [in]  - Matrix (NUMIP x 2) of x-y-coordinates of the vertices cooresponding
01647                             to indices ipindx: \n
01648                             ipcoords(i,0) x-coordinate of node i, \n
01649                             ipcoords(i,1) y-coordinate of node i.
01650   \param  pindx     [in]  - Vector of NUMP indices of ALL nodes in a subdomain, including
01651                             the shared nodes.
01652   \param  pcoords   [in]  - Matrix (NUMP x 2) of x-y-coordinates of the vertices cooresponding
01653                             to indices pindx: \n
01654                             pcoords(i,0) x-coordinate of node i, \n
01655                             pcoords(i,1) y-coordinate of node i.
01656   \param  t         [in]  - Matrix (ELE x 3) of indices of the vertices in a triangle: \n
01657                             t(i,j) index of the j-th vertex in triangle i, where i = 1, ..., ELE
01658   \param  y         [in]  - Reference-counting pointer to the Epetra_MultiVector at which the nonlinear
01659                             term is evaluated.
01660   \param  s         [in]  - Reference-counting pointer to the Epetra_MultiVector that is multiplied by the
01661                             Hessian of the nonlinear form.
01662   \param  lambda    [in]  - Reference-counting pointer to the Epetra_MultiVector of Lagrange Multipliers.
01663   \param  hessvec   [out] - Reference-counting pointer to the Epetra_FECrsMatrix containing Hessian of
01664                             the nonlinear form times vector product.
01665   \return 0                 if successful.
01666 
01667   \par Detailed Description:
01668 
01669   Assembles the nonlinear term \e hessvec, represented by
01670      \f[
01671        \{N''(y,\lambda)s\}_{j} = \langle g''(y_h) \lambda s,\phi_j \rangle
01672                                = \int_{\Omega} g''(y_h(x)) \lambda(x) s(x) \phi_j(x) dx,
01673      \f]
01674   where \f$ g''(y_h) \f$ is given in the local function \e g2pfctn, and \f$\{ \phi_j \}_{j = 1}^{m}\f$ is the
01675   piecewise linear nodal basis for the state space.
01676 */
01677 int GLpApp::nonlinhessvec(const Epetra_Comm & Comm,
01678                   const Epetra_IntSerialDenseVector & ipindx,
01679                   const Epetra_SerialDenseMatrix & ipcoords,
01680                   const Epetra_IntSerialDenseVector & pindx,
01681                   const Epetra_SerialDenseMatrix & pcoords,
01682                   const Epetra_IntSerialDenseMatrix & t,
01683                   const Teuchos::RefCountPtr<const Epetra_MultiVector> & y,
01684                   const Teuchos::RefCountPtr<const Epetra_MultiVector> & s,
01685                   const Teuchos::RefCountPtr<const Epetra_MultiVector> & lambda,
01686                   Teuchos::RefCountPtr<Epetra_FEVector> & hessvec)
01687 {
01688 
01689   int myPID = Comm.MyPID();
01690   int numProcs = Comm.NumProc();
01691 
01692   int numLocNodes     = pindx.M();
01693   int numMyLocNodes   = ipindx.M();
01694   int numLocElems     = t.M();
01695   int numNodesPerElem = 3;
01696 
01697   int indexBase = 1;
01698 
01699   Epetra_Map standardmap(-1, numMyLocNodes, (int*)ipindx.A(), indexBase, Comm);
01700   Epetra_Map overlapmap(-1, numLocNodes, (int*)pindx.A(), indexBase, Comm);
01701 
01702   hessvec = rcp(new Epetra_FEVector(standardmap,1));
01703 
01704   int* nodes = new int[numNodesPerElem];
01705   int i=0, j=0, err=0;
01706   
01707   // get quadrature nodes and weights
01708   Epetra_SerialDenseMatrix Nodes;
01709   Epetra_SerialDenseVector Weights;
01710   quadrature(2,3,Nodes,Weights);
01711   int numQuadPts = Nodes.M();
01712 
01713   // Evaluate nodal basis functions and their derivatives at quadrature points
01714   // N(i,j) = value of the j-th basis function at quadrature node i.
01715   Epetra_SerialDenseMatrix N;
01716   N.Shape(numQuadPts,3);
01717   for (int i=0; i<numQuadPts; i++) {
01718     N(i,0) = 1.0 - Nodes(i,0) - Nodes(i,1);
01719     N(i,1) = Nodes(i,0);
01720     N(i,2) = Nodes(i,1);
01721   }
01722 
01723   // Declare quantities needed for the call to the local assembly routine.
01724   Epetra_IntSerialDenseVector epetra_nodes(View, nodes, numNodesPerElem);
01725   Epetra_SerialDenseMatrix vertices(numNodesPerElem, pcoords.N());
01726 
01727   Epetra_SerialDenseVector ly;        // local entries of y
01728   Epetra_SerialDenseVector Nly;       // N*ly
01729   Epetra_SerialDenseVector ls;        // local entries of s
01730   Epetra_SerialDenseVector Nls;       // N*ls
01731   Epetra_SerialDenseVector llambda;   // local entries of lambda
01732   Epetra_SerialDenseVector Nllambda;  // N*llambda
01733   Epetra_SerialDenseVector lgfctn;    // gfctn(Nly)
01734   Epetra_SerialDenseVector lgfctnNi;  // lgfctn.*N(:,i)
01735   Epetra_SerialDenseVector lgfctnNls; // lgfctn.*N(:,i).*llambda.*ls
01736   Epetra_SerialDenseVector lhessvec;  // local contribution
01737   // Size and init to zero.
01738   ly.Size(numNodesPerElem);
01739   Nly.Size(numQuadPts);
01740   ls.Size(numNodesPerElem);
01741   Nls.Size(numQuadPts);
01742   llambda.Size(numNodesPerElem);
01743   Nllambda.Size(numQuadPts);
01744   lgfctn.Size(numQuadPts);
01745   lgfctnNi.Size(numQuadPts);
01746   lgfctnNls.Size(numQuadPts);
01747   lhessvec.Size(numNodesPerElem);
01748   
01749   Epetra_SerialDenseMatrix B(2,2);
01750   double adB;
01751   
01752   for(i=0; i<numLocElems; i++) {
01753 
01754     nodes[0] = t(i,0); nodes[1] = t(i,1); nodes[2] = t(i,2);
01755     for (j=0; j<numNodesPerElem; j++) {
01756       vertices(j,0) = pcoords(overlapmap.LID(nodes[j]), 0);
01757       vertices(j,1) = pcoords(overlapmap.LID(nodes[j]), 1);
01758     }
01759 
01760     // Construct affine transformation matrix.
01761     for(int i=0; i<2; i++) {
01762       B(i,0) = vertices(1,i)-vertices(0,i);
01763       B(i,1) = vertices(2,i)-vertices(0,i);
01764     }
01765     adB  = abs(determinant(B));
01766 
01767     // Construct local (to each processor) element view of y, s, l.
01768     for (j=0; j<numNodesPerElem; j++) {
01769       ly(j) = (*((*y)(0)))[overlapmap.LID(nodes[j])];
01770       ls(j) = (*((*s)(0)))[overlapmap.LID(nodes[j])];
01771       llambda(j) = (*((*lambda)(0)))[overlapmap.LID(nodes[j])];
01772     }
01773 
01774     Nly.Multiply('N', 'N', 1.0, N, ly, 0.0);            //N*ly
01775     Nls.Multiply('N', 'N', 1.0, N, ls, 0.0);            //N*ls
01776     Nllambda.Multiply('N', 'N', 1.0, N, llambda, 0.0);  //N*llambda
01777     g2pfctn(Nly, lgfctn);
01778     
01779     for (int i=0; i<numNodesPerElem; i++) {
01780       compproduct(lgfctnNi, lgfctn.A(), N[i]);
01781       compproduct(lgfctnNls, lgfctnNi.A(), Nls.A(), Nllambda.A());  // g2pfctn(Nly).*Nls.*Nllambda.*N(:,i)
01782       lhessvec(i) = adB*lgfctnNls.Dot(Weights);
01783     }
01784 
01785     err = hessvec->SumIntoGlobalValues(epetra_nodes, lhessvec);
01786     if (err<0) return(err);
01787   }
01788 
01789   // Call global assemble.
01790 
01791   err = hessvec->GlobalAssemble();
01792   if (err<0) return(err);
01793 
01794   delete [] nodes;
01795 
01796   return(0);
01797 }
01798 
01799 /*  \brief Componentwise evaluation of the second derivative of the nonlinear reaction term.
01800   \param  v   [in]  - Vector at which the second derivative is evaluated.
01801   \param  gv  [out] - Vector value.
01802 */
01803 void GLpApp::g2pfctn(const Epetra_SerialDenseVector & v, Epetra_SerialDenseVector & gv) {
01804   for (int i=0; i<v.M(); i++) {
01805     gv(i) = 6.0*v(i);
01806   }  
01807 }
01808 
01809 /*  \brief Performs finite-element assembly of the Jacobian of the nonlinear form.
01810 
01811   \param  Comm      [in]  - The Epetra (MPI) communicator.
01812   \param  ipindx    [in]  - Vector of NUMIP indices of nodes that are `unique' to a subdomain
01813                             (i.e. owned by the corresponding processor).
01814   \param  ipcoords  [in]  - Matrix (NUMIP x 2) of x-y-coordinates of the vertices cooresponding
01815                             to indices ipindx: \n
01816                             ipcoords(i,0) x-coordinate of node i, \n
01817                             ipcoords(i,1) y-coordinate of node i.
01818   \param  pindx     [in]  - Vector of NUMP indices of ALL nodes in a subdomain, including
01819                             the shared nodes.
01820   \param  pcoords   [in]  - Matrix (NUMP x 2) of x-y-coordinates of the vertices cooresponding
01821                             to indices pindx: \n
01822                             pcoords(i,0) x-coordinate of node i, \n
01823                             pcoords(i,1) y-coordinate of node i.
01824   \param  t         [in]  - Matrix (ELE x 3) of indices of the vertices in a triangle: \n
01825                             t(i,j) index of the j-th vertex in triangle i, where i = 1, ..., ELE
01826   \param  y         [in]  - Reference-counting pointer to the Epetra_MultiVector at which the nonlinear
01827                             term is evaluated.
01828   \param  Gp        [out] - Reference-counting pointer to the Epetra_FECrsMatrix containing the Jacobian
01829                             of the nonlinear form.
01830   \return 0                 if successful.
01831 
01832   \par Detailed Description:
01833 
01834   Assembles the nonlinear term \e Gp, represented by
01835      \f[
01836        \{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,
01837      \f]
01838   where \f$ g'(y_h) \f$ is given in the local function \e gpfctn, and \f$\{ \phi_j \}_{j = 1}^{m}\f$ is the
01839   piecewise linear nodal basis for the state space.
01840 */
01841 int GLpApp::nonlinjac(const Epetra_Comm & Comm,
01842               const Epetra_IntSerialDenseVector & ipindx,
01843               const Epetra_SerialDenseMatrix & ipcoords,
01844               const Epetra_IntSerialDenseVector & pindx,
01845               const Epetra_SerialDenseMatrix & pcoords,
01846               const Epetra_IntSerialDenseMatrix & t,
01847               const Teuchos::RefCountPtr<const Epetra_MultiVector> & y,
01848               Teuchos::RefCountPtr<Epetra_FECrsMatrix> & Gp)
01849 {
01850 
01851   int myPID = Comm.MyPID();
01852   int numProcs = Comm.NumProc();
01853 
01854   int numLocNodes     = pindx.M();
01855   int numMyLocNodes   = ipindx.M();
01856   int numLocElems     = t.M();
01857   int numNodesPerElem = 3;
01858 
01859   int indexBase = 1;
01860 
01861   Epetra_Map standardmap(-1, numMyLocNodes, (int*)ipindx.A(), indexBase, Comm);
01862   Epetra_Map overlapmap(-1, numLocNodes, (int*)pindx.A(), indexBase, Comm);
01863 
01864   int format = Epetra_FECrsMatrix::COLUMN_MAJOR;
01865   Gp = rcp(new Epetra_FECrsMatrix(Copy, standardmap, 0));
01866 
01867   int* nodes = new int[numNodesPerElem];
01868   int i=0, j=0, err=0;
01869   
01870   // get quadrature nodes and weights
01871   Epetra_SerialDenseMatrix Nodes;
01872   Epetra_SerialDenseVector Weights;
01873   quadrature(2,3,Nodes,Weights);
01874   int numQuadPts = Nodes.M();
01875 
01876   // Evaluate nodal basis functions and their derivatives at quadrature points
01877   // N(i,j) = value of the j-th basis function at quadrature node i.
01878   Epetra_SerialDenseMatrix N;
01879   N.Shape(numQuadPts,3);
01880   for (int i=0; i<numQuadPts; i++) {
01881     N(i,0) = 1.0 - Nodes(i,0) - Nodes(i,1);
01882     N(i,1) = Nodes(i,0);
01883     N(i,2) = Nodes(i,1);
01884   }
01885   
01886   // Declare quantities needed for the call to the local assembly routine.
01887   Epetra_IntSerialDenseVector epetra_nodes(View, nodes, numNodesPerElem);
01888   Epetra_SerialDenseMatrix vertices(numNodesPerElem, pcoords.N());
01889 
01890   Epetra_SerialDenseVector ly;          // local entries of y
01891   Epetra_SerialDenseVector Nly;         // N*ly
01892   Epetra_SerialDenseVector lgfctn;      // gfctn(Nly)
01893   Epetra_SerialDenseVector lgfctnNiNj;  // lgfctn.*N(:,i).*N(:,j)
01894   Epetra_SerialDenseMatrix lGp;         // local contribution
01895   // Size and init to zero.
01896   ly.Size(numNodesPerElem);
01897   Nly.Size(numQuadPts);
01898   lgfctn.Size(numQuadPts);
01899   lgfctnNiNj.Size(numQuadPts);
01900   lGp.Shape(numNodesPerElem, numNodesPerElem);
01901   
01902   Epetra_SerialDenseMatrix B(2,2);
01903   double adB;
01904   
01905   for(i=0; i<numLocElems; i++) {
01906 
01907     nodes[0] = t(i,0); nodes[1] = t(i,1); nodes[2] = t(i,2);
01908     for (j=0; j<numNodesPerElem; j++) {
01909       vertices(j,0) = pcoords(overlapmap.LID(nodes[j]), 0);
01910       vertices(j,1) = pcoords(overlapmap.LID(nodes[j]), 1);
01911     }
01912 
01913     // Construct affine transformation matrix.
01914     for(int i=0; i<2; i++) {
01915       B(i,0) = vertices(1,i)-vertices(0,i);
01916       B(i,1) = vertices(2,i)-vertices(0,i);
01917     }
01918     adB  = abs(determinant(B));
01919 
01920     // Construct local (to each processor) element view of y. 
01921     for (j=0; j<numNodesPerElem; j++) {
01922       ly(j) = (*((*y)(0)))[overlapmap.LID(nodes[j])];
01923     }
01924 
01925     Nly.Multiply('N', 'N', 1.0, N, ly, 0.0);
01926     gpfctn(Nly, lgfctn);
01927     
01928     for (int i=0; i<numNodesPerElem; i++) {
01929       for (int j=0; j<numNodesPerElem; j++) {
01930         compproduct(lgfctnNiNj, lgfctn.A(), N[i], N[j]);
01931         lGp(i,j) = adB*lgfctnNiNj.Dot(Weights);
01932       }
01933     }
01934     
01935     err = Gp->InsertGlobalValues(epetra_nodes, lGp, format);
01936     if (err<0) return(err);
01937   }
01938 
01939   // Call global assemble.
01940 
01941   err = Gp->GlobalAssemble();
01942   if (err<0) return(err);
01943 
01944   delete [] nodes;
01945 
01946   return(0);
01947 }
01948 
01949 /*  \brief Componentwise evaluation of the first derivative of the nonlinear reaction term.
01950   \param  v   [in]  - Vector at which the first derivative is evaluated.
01951   \param  gv  [out] - Vector value.
01952 */
01953 void GLpApp::gpfctn(const Epetra_SerialDenseVector & v, Epetra_SerialDenseVector & gv) {
01954   for (int i=0; i<v.M(); i++) {
01955     gv(i) = 3.0*pow(v(i),2)-1.0;
01956   }  
01957 }
01958 
01959 /*  \brief Performs finite-element assembly of the nonlinear reaction term.
01960 
01961   \param  Comm      [in]  - The Epetra (MPI) communicator.
01962   \param  ipindx    [in]  - Vector of NUMIP indices of nodes that are `unique' to a subdomain
01963                             (i.e. owned by the corresponding processor).
01964   \param  ipcoords  [in]  - Matrix (NUMIP x 2) of x-y-coordinates of the vertices cooresponding
01965                             to indices ipindx: \n
01966                             ipcoords(i,0) x-coordinate of node i, \n
01967                             ipcoords(i,1) y-coordinate of node i.
01968   \param  pindx     [in]  - Vector of NUMP indices of ALL nodes in a subdomain, including
01969                             the shared nodes.
01970   \param  pcoords   [in]  - Matrix (NUMP x 2) of x-y-coordinates of the vertices cooresponding
01971                             to indices pindx: \n
01972                             pcoords(i,0) x-coordinate of node i, \n
01973                             pcoords(i,1) y-coordinate of node i.
01974   \param  t         [in]  - Matrix (ELE x 3) of indices of the vertices in a triangle: \n
01975                             t(i,j) index of the j-th vertex in triangle i, where i = 1, ..., ELE
01976   \param  y         [out] - Reference-counting pointer to the Epetra_MultiVector at which the nonlinear
01977                             form is evaluated.
01978   \param  g         [out] - Reference-counting pointer to the Epetra_FEVector containing the value
01979                             of the nonlinear form.
01980   \return 0                 if successful.
01981 
01982   \par Detailed Description:
01983 
01984   Assembles the nonlinear term \e g, represented by
01985      \f[
01986        \{N(y)\}_{j} = \langle g(y_h),\phi_j \rangle =  \int_{\Omega} g(y_h(x)) \phi_j(x) dx,
01987      \f]
01988   where \f$ g(y_h) \f$ is given in the local function \e gfctn, and \f$\{ \phi_j \}_{j = 1}^{m}\f$ is the
01989   piecewise linear nodal basis for the state space.
01990 */
01991 int GLpApp::nonlinvec(const Epetra_Comm & Comm,
01992               const Epetra_IntSerialDenseVector & ipindx,
01993               const Epetra_SerialDenseMatrix & ipcoords,
01994               const Epetra_IntSerialDenseVector & pindx,
01995               const Epetra_SerialDenseMatrix & pcoords,
01996               const Epetra_IntSerialDenseMatrix & t,
01997               const Teuchos::RefCountPtr<const Epetra_MultiVector> & y,
01998               Teuchos::RefCountPtr<Epetra_FEVector> & g)
01999 {
02000 
02001   int myPID = Comm.MyPID();
02002   int numProcs = Comm.NumProc();
02003 
02004   int numLocNodes     = pindx.M();
02005   int numMyLocNodes   = ipindx.M();
02006   int numLocElems     = t.M();
02007   int numNodesPerElem = 3;
02008 
02009   int indexBase = 1;
02010 
02011   Epetra_Map standardmap(-1, numMyLocNodes, (int*)ipindx.A(), indexBase, Comm);
02012   Epetra_Map overlapmap(-1, numLocNodes, (int*)pindx.A(), indexBase, Comm);
02013 
02014   g = rcp(new Epetra_FEVector(standardmap,1));
02015 
02016   int* nodes = new int[numNodesPerElem];
02017   int i=0, j=0, err=0;
02018   
02019   // get quadrature nodes and weights
02020   Epetra_SerialDenseMatrix Nodes;
02021   Epetra_SerialDenseVector Weights;
02022   quadrature(2,3,Nodes,Weights);
02023   int numQuadPts = Nodes.M();
02024 
02025   // Evaluate nodal basis functions and their derivatives at quadrature points
02026   // N(i,j) = value of the j-th basis function at quadrature node i.
02027   Epetra_SerialDenseMatrix N;
02028   N.Shape(numQuadPts,3);
02029   for (int i=0; i<numQuadPts; i++) {
02030     N(i,0) = 1.0 - Nodes(i,0) - Nodes(i,1);
02031     N(i,1) = Nodes(i,0);
02032     N(i,2) = Nodes(i,1);
02033   }
02034 
02035   // Declare quantities needed for the call to the local assembly routine.
02036   Epetra_IntSerialDenseVector epetra_nodes(View, nodes, numNodesPerElem);
02037   Epetra_SerialDenseMatrix vertices(numNodesPerElem, pcoords.N());
02038 
02039   Epetra_SerialDenseVector ly;        // local entries of y
02040   Epetra_SerialDenseVector Nly;       // N*ly
02041   Epetra_SerialDenseVector lgfctn;    // gfctn(Nly)
02042   Epetra_SerialDenseVector lgfctnNi;  // lgfctn.*N(:,i)
02043   Epetra_SerialDenseVector lg;        // local contribution
02044   // Size and init to zero.
02045   ly.Size(numNodesPerElem);
02046   Nly.Size(numQuadPts);
02047   lgfctn.Size(numQuadPts);
02048   lgfctnNi.Size(numQuadPts);
02049   lg.Size(numNodesPerElem);
02050   
02051   Epetra_SerialDenseMatrix B(2,2);
02052   double adB;
02053   
02054   for(i=0; i<numLocElems; i++) {
02055 
02056     nodes[0] = t(i,0); nodes[1] = t(i,1); nodes[2] = t(i,2);
02057     for (j=0; j<numNodesPerElem; j++) {
02058       vertices(j,0) = pcoords(overlapmap.LID(nodes[j]), 0);
02059       vertices(j,1) = pcoords(overlapmap.LID(nodes[j]), 1);
02060     }
02061 
02062     // Construct affine transformation matrix.
02063     for(int i=0; i<2; i++) {
02064       B(i,0) = vertices(1,i)-vertices(0,i);
02065       B(i,1) = vertices(2,i)-vertices(0,i);
02066     }
02067     adB  = abs(determinant(B));
02068 
02069     // Construct local (to each processor) element view of y. 
02070     for (j=0; j<numNodesPerElem; j++) {
02071       ly(j) = (*((*y)(0)))[overlapmap.LID(nodes[j])];
02072     }
02073 
02074     Nly.Multiply('N', 'N', 1.0, N, ly, 0.0);
02075     gfctn(Nly, lgfctn);
02076     
02077     for (int i=0; i<numNodesPerElem; i++) {
02078       compproduct(lgfctnNi, lgfctn.A(), N[i]);
02079       lg(i) = adB*lgfctnNi.Dot(Weights);
02080     }
02081     
02082     err = g->SumIntoGlobalValues(epetra_nodes, lg);
02083     if (err<0) return(err);
02084   }
02085 
02086   // Call global assemble.
02087 
02088   err = g->GlobalAssemble();
02089   if (err<0) return(err);
02090 
02091   delete [] nodes;
02092 
02093   return(0);
02094 }
02095 
02096 
02097 /*  \brief Componentwise evaluation of the nonlinear reaction term.
02098   \param  v   [in]  - Vector at which the nonlinear function is evaluated.
02099   \param  gv  [out] - Vector value.
02100 */
02101 void GLpApp::gfctn(const Epetra_SerialDenseVector & v, Epetra_SerialDenseVector & gv) {
02102   for (int i=0; i<v.M(); i++) {
02103     gv(i) = pow(v(i),3)-v(i);
02104   }  
02105 }
02106 
02107 /* ======== ================ *
02108  * function CrsMatrix2MATLAB *
02109  * ======== ================ *
02110  *
02111  * Print out a CrsMatrix in a MATLAB format. Each processor prints out
02112  * its part, starting from proc 0 to proc NumProc-1. The first line of
02113  * each processor's output states the number of local rows and of
02114  * local nonzero elements.
02115  *
02116  *
02117  * Return code:        true if matrix has been printed out
02118  * -----------         false otherwise
02119  *
02120  * Parameters:
02121  * ----------
02122  *
02123  * - Epetra_CrsMatrix  reference to the distributed CrsMatrix to 
02124  *                     print out
02125  * - ostream &         reference to output stream 
02126  */
02127 
02128 bool GLpApp::CrsMatrix2MATLAB(const Epetra_CrsMatrix & A, ostream & outfile) 
02129 {
02130 
02131   int MyPID = A.Comm().MyPID(); 
02132   int NumProc = A.Comm().NumProc();
02133 
02134   // work only on transformed matrices;
02135   if( A.IndicesAreLocal() == false ) {
02136     if( MyPID == 0 ) { 
02137       cerr << "ERROR in "<< __FILE__ << ", line " << __LINE__ << endl;
02138       cerr << "Function CrsMatrix2MATLAB accepts\n";
02139       cerr << "transformed matrices ONLY. Please call A.TransformToLoca()\n";
02140       cerr << "on your matrix A to that purpose.\n";
02141       cerr << "Now returning...\n";
02142     }
02143     return false;
02144   }
02145 
02146   int NumMyRows = A.NumMyRows(); // number of rows on this process
02147   int NumNzRow;   // number of nonzero elements for each row
02148   int NumEntries; // number of extracted elements for each row
02149   int NumGlobalRows; // global dimensio of the problem
02150   int GlobalRow;  // row in global ordering
02151   int NumGlobalNonzeros; // global number of nonzero elements
02152 
02153   NumGlobalRows = A.NumGlobalRows();
02154   NumGlobalNonzeros = A.NumGlobalNonzeros();
02155 
02156   // print out on cout if no filename is provided
02157 
02158   int IndexBase = A.IndexBase(); // MATLAB starts from 0
02159   if( IndexBase == 0 )
02160     IndexBase = 1;
02161   else if ( IndexBase == 1)
02162     IndexBase = 0;
02163 
02164   // write on file the dimension of the matrix
02165 
02166   if( MyPID==0 ) {
02167     outfile << "A = spalloc(";
02168     outfile << NumGlobalRows << ',' << NumGlobalRows;
02169     outfile << ',' << NumGlobalNonzeros << ");\n";
02170   }
02171 
02172   for( int Proc=0 ; Proc<NumProc ; ++Proc ) {
02173     A.Comm().Barrier();
02174     if( MyPID == Proc ) {
02175 
02176       outfile << "\n\n% On proc " << Proc << ": ";
02177       outfile << NumMyRows << " rows and ";
02178       outfile << A.NumMyNonzeros() << " nonzeros\n";
02179 
02180       // cycle over all local rows to find out nonzero elements
02181       for( int MyRow=0 ; MyRow<NumMyRows ; ++MyRow ) {
02182 
02183         GlobalRow = A.GRID(MyRow);
02184 
02185         NumNzRow = A.NumMyEntries(MyRow);
02186         double *Values = new double[NumNzRow];
02187         int *Indices = new int[NumNzRow];
02188 
02189         A.ExtractMyRowCopy(MyRow, NumNzRow, 
02190                            NumEntries, Values, Indices);
02191         // print out the elements with MATLAB syntax
02192         for( int j=0 ; j<NumEntries ; ++j ) {
02193           outfile << "A(" << GlobalRow  + IndexBase 
02194                   << "," << A.GCID(Indices[j]) + IndexBase
02195                   << ") = " << Values[j] << ";\n";
02196         }
02197 
02198         delete Values;
02199         delete Indices;
02200       }
02201       
02202     }
02203     A.Comm().Barrier();
02204   }
02205 
02206   return true;
02207 
02208 }
02209 
02210 
02211 /* ======== ============= *
02212  * function Vector2MATLAB *
02213  * ======== ============= *
02214  *
02215  * Print out a Epetra_Vector in a MATLAB format. Each processor prints out
02216  * its part, starting from proc 0 to proc NumProc-1. The first line of
02217  * each processor's output states the number of local rows and of
02218  * local nonzero elements.
02219  *
02220  * Return code:        true if vector has been printed out
02221  * -----------         false otherwise
02222  *
02223  * Parameters:
02224  * ----------
02225  *
02226  * - Epetra_Vector     reference to vector
02227  * - ostream &         reference to output stream 
02228  */
02229 
02230 bool GLpApp::Vector2MATLAB( const Epetra_Vector & v, ostream & outfile)
02231 {
02232   
02233   int MyPID = v.Comm().MyPID(); 
02234   int NumProc = v.Comm().NumProc();
02235   int MyLength = v.MyLength();
02236   int GlobalLength = v.GlobalLength();
02237   
02238   // print out on cout if no filename is provided
02239 
02240   // write on file the dimension of the matrix
02241 
02242   if( MyPID == 0 ) outfile << "v = zeros(" << GlobalLength << ",1)\n";
02243 
02244   int NumMyElements = v.Map().NumMyElements();
02245   // get update list
02246   int * MyGlobalElements = v.Map().MyGlobalElements( );
02247   
02248   int Row;
02249 
02250   int IndexBase = v.Map().IndexBase(); // MATLAB starts from 0
02251   if( IndexBase == 0 )
02252     IndexBase = 1;
02253   else if ( IndexBase == 1)
02254     IndexBase = 0;
02255 
02256   for( int Proc=0 ; Proc<NumProc ; ++Proc ) {
02257     v.Comm().Barrier();
02258     if( MyPID == Proc ) {
02259 
02260       outfile << "% On proc " << Proc << ": ";
02261       outfile << MyLength << " rows of ";
02262       outfile << GlobalLength << " elements\n";
02263 
02264       for( Row=0 ; Row<MyLength ; ++Row ) {
02265         outfile << "v(" << MyGlobalElements[Row] + IndexBase
02266              << ") = " << v[Row] << ";\n";
02267       }
02268       
02269     }
02270       
02271     v.Comm().Barrier();
02272   }
02273 
02274   return true;
02275 
02276 } /* Vector2MATLAB */
02277 
02278 
02279 /* ======== =============== *
02280  * function FEVector2MATLAB *
02281  * ======== =============== *
02282  *
02283  * Print out a Epetra_Vector in a MATLAB format. Each processor prints out
02284  * its part, starting from proc 0 to proc NumProc-1. The first line of
02285  * each processor's output states the number of local rows and of
02286  * local nonzero elements.
02287  *
02288  * Return code:        true if vector has been printed out
02289  * -----------         false otherwise
02290  *
02291  * Parameters:
02292  * ----------
02293  *
02294  * - Epetra_FEVector   reference to FE vector
02295  * - ostream &         reference to output stream 
02296  */
02297 
02298 bool GLpApp::FEVector2MATLAB( const Epetra_FEVector & v, ostream & outfile)
02299 {
02300   
02301   int MyPID = v.Comm().MyPID(); 
02302   int NumProc = v.Comm().NumProc();
02303   int MyLength = v.MyLength();
02304   int GlobalLength = v.GlobalLength();
02305   
02306   // print out on cout if no filename is provided
02307 
02308   // write on file the dimension of the matrix
02309 
02310   if( MyPID == 0 ) outfile << "v = zeros(" << GlobalLength << ",1)\n";
02311 
02312   int NumMyElements = v.Map().NumMyElements();
02313   // get update list
02314   int * MyGlobalElements = v.Map().MyGlobalElements( );
02315   
02316   int Row;
02317 
02318   int IndexBase = v.Map().IndexBase(); // MATLAB starts from 0
02319   if( IndexBase == 0 )
02320     IndexBase = 1;
02321   else if ( IndexBase == 1)
02322     IndexBase = 0;
02323 
02324   for( int Proc=0 ; Proc<NumProc ; ++Proc ) {
02325     v.Comm().Barrier();
02326     if( MyPID == Proc ) {
02327 
02328       outfile << "% On proc " << Proc << ": ";
02329       outfile << MyLength << " rows of ";
02330       outfile << GlobalLength << " elements\n";
02331 
02332       for( Row=0 ; Row<MyLength ; ++Row ) {
02333         outfile << "v(" << MyGlobalElements[Row] + IndexBase
02334              << ") = " << v[0][Row] << ";\n";
02335       }
02336       
02337     }
02338       
02339     v.Comm().Barrier();
02340   }
02341 
02342   return true;
02343 
02344 } /* FEVector2MATLAB */
02345 
02346 
02347 /*  \brief  Returns the nodes and weights for the integration \n
02348              on the interval [0,1] (dim = 1) \n
02349              on the triangle with vertices (0,0), (1,0), (0,1) (if dim = 2) \n
02350              on the tetrahedron with vertices (0,0,0), (1,0,0), (0,1,0), (0,0,1) (if dim = 3).
02351              
02352   \param  dim     [in]   - spatial dimension (dim = 1, 2)
02353   \param  order   [in]   - required degree of polynomials that integrate exactly
02354   \param  nodes   [out]  - Matrix in which the i-th row of nodes gives coordinates of the i-th quadrature node
02355   \param  weights [out]  - quadrature weights
02356 
02357   \return 0                if successful
02358 */
02359 int GLpApp::quadrature(const int dim, const int order,
02360                Epetra_SerialDenseMatrix & nodes,
02361                Epetra_SerialDenseVector & weights)
02362 {
02363   
02364   if (dim == 1) {
02365 
02366     // Gauss quadrature nodes and weights on the interval [0,1]
02367 
02368     if (order == 1) {
02369       nodes.Shape(1,1);
02370       nodes(0,0) = 0.5;
02371       weights.Size(1);
02372       weights(0) = 1.0;
02373     }
02374     else if (order == 2) {
02375       nodes.Shape(2,1);
02376       nodes(0,0) = (1.0-1.0/sqrt(3.0))/2.0;
02377       nodes(1,0) = (1.0+1.0/sqrt(3.0))/2.0;
02378       weights.Size(2);
02379       weights(0) = 0.5;
02380       weights(1) = 0.5;
02381     }
02382     else if (order == 3) {
02383       nodes.Shape(3,1);
02384       nodes(0,0) = (1.0-sqrt(3.0/5.0))/2.0;
02385       nodes(1,0) = 0.5;
02386       nodes(2,0) = (1.0+sqrt(3.0/5.0))/2.0;
02387       weights.Size(3);
02388       weights(0) = 5.0/18.0;
02389       weights(1) = 4.0/9.0;
02390       weights(2) = 5.0/18.0;
02391     }
02392     else {
02393       cout << "Quadrature for dim = " << dim << " and order = ";
02394       cout << order << " not available.\n";
02395       exit(-1);
02396     }
02397 
02398   }
02399   else if (dim == 2) {
02400     
02401     // Quadrature nodes and weights on the unit simplex with
02402     // vertices (0,0), (1,0), and (0,1).
02403 
02404     if (order == 1) {
02405       nodes.Shape(1,2);
02406       nodes(0,0) = 1.0/3.0; nodes (0,1) = 1.0/3.0;
02407       weights.Size(1);
02408       weights(0) = 0.5;
02409     }
02410     else if (order == 2) {
02411       nodes.Shape(3,2);
02412       nodes(0,0) = 1.0/6.0; nodes (0,1) = 1.0/6.0;
02413       nodes(1,0) = 2.0/3.0; nodes (1,1) = 1.0/6.0;
02414       nodes(2,0) = 1.0/6.0; nodes (2,1) = 2.0/3.0;
02415       weights.Size(3);
02416       weights(0) = 1.0/6.0;
02417       weights(1) = 1.0/6.0;
02418       weights(2) = 1.0/6.0;
02419     }
02420     else if (order == 3) {
02421       nodes.Shape(4,2);
02422       nodes(0,0) = 1.0/3.0; nodes (0,1) = 1.0/3.0;
02423       nodes(1,0) = 3.0/5.0; nodes (1,1) = 1.0/5.0;
02424       nodes(2,0) = 1.0/5.0; nodes (2,1) = 3.0/5.0;
02425       nodes(3,0) = 1.0/5.0; nodes (3,1) = 1.0/5.0;
02426       weights.Size(4);
02427       weights(0) = -9.0/32.0;
02428       weights(1) = 25.0/96.0;
02429       weights(2) = 25.0/96.0;
02430       weights(3) = 25.0/96.0;
02431     }
02432     else {
02433       cout << "Quadrature for dim = " << dim << " and order = ";
02434       cout << order << " not available.\n";
02435       exit(-1);
02436     }
02437 
02438   }
02439   else {
02440     cout << "Quadrature for dim = " << dim << " not available.\n";
02441     exit(-1);
02442   }
02443 
02444   return(0);
02445 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines