MueLu_Q2Q1Q2CoarseGridFactory_def.hpp

Go to the documentation of this file.

// @HEADER
//
// ***********************************************************************
//
//        MueLu: A package for multigrid based preconditioning
//                  Copyright 2012 Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact
//                    Jonathan Hu       (jhu@sandia.gov)
//                    Andrey Prokopenko (aprokop@sandia.gov)
//                    Ray Tuminaro      (rstumin@sandia.gov)
//
// ***********************************************************************
//
// @HEADER
#ifndef MUELU_Q2Q1Q2COARSEGRIDFACTORY_DEF_HPP
#define MUELU_Q2Q1Q2COARSEGRIDFACTORY_DEF_HPP

#include <iostream>
#include <cmath>

#include <Teuchos_SerialDenseMatrix.hpp>

#include <Xpetra_Map.hpp>
#include <Xpetra_MapFactory.hpp>
#include <Xpetra_Matrix.hpp>
#include <Xpetra_MultiVector.hpp>
#include <Xpetra_MultiVectorFactory.hpp>
#include <Xpetra_VectorFactory.hpp>
#include <Xpetra_CrsMatrixWrap.hpp>
#include <Xpetra_MultiVector.hpp>
#include <Xpetra_MultiVectorFactory.hpp>

#include "MueLu_Q2Q1Q2CoarseGridFactory_decl.hpp"
#include <MueLu_Level.hpp>
#include <MueLu_Monitor.hpp>
#include <MueLu_Utilities.hpp>
#include <MueLu_FactoryManagerBase.hpp>
#include "MueLu_SingleLevelFactoryBase.hpp"

namespace MueLu {

  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  Q2Q1Q2CoarseGridFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Q2Q1Q2CoarseGridFactory(){
    GetOStream(Runtime1) << "I constructed a Q2Q1Q2CoarseGridFactory object... Nothing else to do here." << std::endl;
  }

  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  Q2Q1Q2CoarseGridFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::~Q2Q1Q2CoarseGridFactory(){
    // Should be empty. All destruction should be handled by Level-based get stuff and RCP
  }

  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  void Q2Q1Q2CoarseGridFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level &fineLevel, Level &coarseLevel) const {

    Input(fineLevel, "VElementList");
    Input(fineLevel, "PElementList");
    Input(fineLevel, "MElementList");

    Input(coarseLevel, "VElementList");
    Input(coarseLevel, "PElementList");
    Input(coarseLevel, "MElementList");


    //currentLevel.DeclareInput(varName_,factory_,this);
  }

  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  void Q2Q1Q2CoarseGridFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level &fineLevel, Level &coarseLevel) const {

    GetOStream(Runtime1) << "Starting 'build' routine...\n";

    // This will create a list of elements on the coarse grid with a
    // predictable structure, as well as modify the fine grid list of
    // elements, if necessary (i.e. if fineLevel.GetLevelID()==0);
    //BuildCoarseGrid(fineLevel,coarseLevel);

    // This will actually build our prolongator P
    return BuildCoarseGrid(fineLevel,coarseLevel);

  }


  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  void Q2Q1Q2CoarseGridFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildCoarseGrid(Level &fineLevel, Level &coarseLevel) const
  {

    GetOStream(Runtime1) << "starting 'BuildCoarseGrid' routine...\n";

    RCP<Teuchos::SerialDenseMatrix<GO,GO> > fineElementPDOFs = Get< RCP<Teuchos::SerialDenseMatrix<GO,GO> > >(fineLevel,"PElementList");

    GO totalFineElements = fineElementPDOFs->numRows();

    // Compute number of coarse grid elements in total:
    GO totalCoarseElements = totalFineElements/4;
    LO nCoarseElements = (int) sqrt(totalCoarseElements);

    // Initialize some counters:
    size_t EdgeCount = (nCoarseElements + 1) * (nCoarseElements + 1);
    size_t CenterCount = EdgeCount + 2 * nCoarseElements * (nCoarseElements + 1);


    // Initialize arrays of the proper size:
    RCP<Teuchos::SerialDenseMatrix<GO,GO> > coarseElementVDOFs = rcp(new Teuchos::SerialDenseMatrix<GO,GO>(totalCoarseElements,18));
    RCP<Teuchos::SerialDenseMatrix<GO,GO> > coarseElementPDOFs = rcp(new Teuchos::SerialDenseMatrix<GO,GO>(totalCoarseElements,4));
    RCP<Teuchos::SerialDenseMatrix<GO,GO> > coarseElementMDOFs = rcp(new Teuchos::SerialDenseMatrix<GO,GO>(totalCoarseElements,9));


    for ( GO coarseElement=0; coarseElement < totalCoarseElements; coarseElement++ )
      {

        // ***************************************************************
        // This is less of a pain in the ass for magnetics, so I'm
        // going to build the magnetics list. The velocity follows
        // by doubling everything (and adding 1 for the y-components)
        // and the pressure follows by copying the magnetics nodes.
        // ***************************************************************

        //if (coarseElement is on the Bottom Edge)
        if (coarseElement < nCoarseElements)
          {
            //Bottom nodes
            (*coarseElementMDOFs)(coarseElement,0) = coarseElement;
            (*coarseElementMDOFs)(coarseElement,1) = coarseElement+1;

            //Bottom edge
            (*coarseElementMDOFs)(coarseElement,4) = EdgeCount++;

          }
        else
          {
            //Bottom Nodes
            (*coarseElementMDOFs)(coarseElement,0) = (*coarseElementMDOFs)(coarseElement-nCoarseElements,3);
            (*coarseElementMDOFs)(coarseElement,1) = (*coarseElementMDOFs)(coarseElement-nCoarseElements,2);

            //Bottom Edge
            (*coarseElementMDOFs)(coarseElement,4) = (*coarseElementMDOFs)(coarseElement-nCoarseElements,6);

          }


        //Right and Top Edges -- must be determined before left edge
        (*coarseElementMDOFs)(coarseElement,5) = EdgeCount++;
        (*coarseElementMDOFs)(coarseElement,6) = EdgeCount++;


        //if (coarseElement is on the Left Edge)
        if (coarseElement % nCoarseElements == 0)
          {
            //Top left node
            (*coarseElementMDOFs)(coarseElement,3) = (*coarseElementMDOFs)(coarseElement,0)+nCoarseElements+1;

            //Left Edge
            (*coarseElementMDOFs)(coarseElement,7) = EdgeCount++;

          }
        else
          {
            //Top left node
            (*coarseElementMDOFs)(coarseElement,3) = (*coarseElementMDOFs)(coarseElement-1,2);

            //Left Edge
            (*coarseElementMDOFs)(coarseElement,7) = (*coarseElementMDOFs)(coarseElement-1,5);
          }

        //Top right node -- Must be the last node to be determined!
        (*coarseElementMDOFs)(coarseElement,2) = (*coarseElementMDOFs)(coarseElement,3)+1;

        //Center Node
        (*coarseElementMDOFs)(coarseElement,8) = CenterCount++;



        // With Magnetics built, Pressure and Velocity follow without effort.
        // First, Velocity:
        (*coarseElementVDOFs)(coarseElement,0) = 2*(*coarseElementMDOFs)(coarseElement,0);
        (*coarseElementVDOFs)(coarseElement,1) = 2*(*coarseElementMDOFs)(coarseElement,0)+1;
        (*coarseElementVDOFs)(coarseElement,2) = 2*(*coarseElementMDOFs)(coarseElement,1);
        (*coarseElementVDOFs)(coarseElement,3) = 2*(*coarseElementMDOFs)(coarseElement,1)+1;
        (*coarseElementVDOFs)(coarseElement,4) = 2*(*coarseElementMDOFs)(coarseElement,2);
        (*coarseElementVDOFs)(coarseElement,5) = 2*(*coarseElementMDOFs)(coarseElement,2)+1;
        (*coarseElementVDOFs)(coarseElement,6) = 2*(*coarseElementMDOFs)(coarseElement,3);
        (*coarseElementVDOFs)(coarseElement,7) = 2*(*coarseElementMDOFs)(coarseElement,3)+1;
        (*coarseElementVDOFs)(coarseElement,8) = 2*(*coarseElementMDOFs)(coarseElement,4);
        (*coarseElementVDOFs)(coarseElement,9) = 2*(*coarseElementMDOFs)(coarseElement,4)+1;
        (*coarseElementVDOFs)(coarseElement,10) = 2*(*coarseElementMDOFs)(coarseElement,5);
        (*coarseElementVDOFs)(coarseElement,11) = 2*(*coarseElementMDOFs)(coarseElement,5)+1;
        (*coarseElementVDOFs)(coarseElement,12) = 2*(*coarseElementMDOFs)(coarseElement,6);
        (*coarseElementVDOFs)(coarseElement,13) = 2*(*coarseElementMDOFs)(coarseElement,6)+1;
        (*coarseElementVDOFs)(coarseElement,14) = 2*(*coarseElementMDOFs)(coarseElement,7);
        (*coarseElementVDOFs)(coarseElement,15) = 2*(*coarseElementMDOFs)(coarseElement,7)+1;
        (*coarseElementVDOFs)(coarseElement,16) = 2*(*coarseElementMDOFs)(coarseElement,8);
        (*coarseElementVDOFs)(coarseElement,17) = 2*(*coarseElementMDOFs)(coarseElement,8)+1;

        // Lastly, Pressure:
        (*coarseElementPDOFs)(coarseElement,0) = (*coarseElementMDOFs)(coarseElement,0);
        (*coarseElementPDOFs)(coarseElement,1) = (*coarseElementMDOFs)(coarseElement,1);
        (*coarseElementPDOFs)(coarseElement,2) = (*coarseElementMDOFs)(coarseElement,2);
        (*coarseElementPDOFs)(coarseElement,3) = (*coarseElementMDOFs)(coarseElement,3);

      }// Loop over elements

    Set(coarseLevel,"VElementList",coarseElementVDOFs);
    Set(coarseLevel,"PElementList",coarseElementPDOFs);
    Set(coarseLevel,"MElementList",coarseElementMDOFs);

    //coarseLevel.Keep("VElementList",coarseLevel.GetFactoryManager()->GetFactory("VElementList").get());
    //coarseLevel.Keep("PElementList",coarseLevel.GetFactoryManager()->GetFactory("PElementList").get());
    //coarseLevel.Keep("MElementList",coarseLevel.GetFactoryManager()->GetFactory("MElementList").get());

  }//BuildCoarseGrid


} // namespace MueLu

#define MUELU_Q2Q1Q2COARSEGRIDFACTORY_SHORT
#endif // MUELU_Q2Q1Q2COARSEGRIDFACTORY_DEF_HPP