TPZRefPattern.cpp

Go to the documentation of this file.

// /**
//  * @file
//  * @brief Contains the implementation of the TPZRefPattern methods.
//  */
// /* Generated by Together */

#include "TPZRefPattern.h"
#include "TPZRefPatternTools.h"
#include "TPZRefPatternDataBase.h"
#include "TPZVTKGeoMesh.h"
// #include "pztrnsform.h"
// #include "pzreal.h"
// #include "pzgmesh.h"
// #include "pzquad.h"
// #include "pzvec.h"
// #include "pzeltype.h"
// #include "tpzpermutation.h"
// #include "pzgeoel.h"
// #include "pzlog.h"

// #include <set>

// #include <fstream>
// #include <sstream>

//iniatilizing static variables
const int TPZRefPattern::fNonInitializedId = -50;
const std::string TPZRefPattern::fNonInitializedName = "noname";
std::map<MElementType, std::list<TPZRefPattern::TPZRefPatternPermute> > TPZRefPattern::fPermutations;

 #ifdef LOG4CXX
 static LoggerPtr logger(Logger::getLogger("pz.mesh.TPZRefPattern"));
 #endif


TPZRefPattern::TPZRefPattern() : fId(fNonInitializedId), fName(fNonInitializedName),
fSideRefPattern(0), fPermutedRefPatterns(0), fNSubEl(0){

}

TPZRefPattern::TPZRefPattern(const TPZRefPattern &cp){
    fSubElSideInfo = cp.fSubElSideInfo;
    fFatherSideInfo = cp.fFatherSideInfo;
    fName = cp.fName;
    fId = cp.fId;
    fSideRefPattern = cp.fSideRefPattern;
    fPermutedRefPatterns = cp.fPermutedRefPatterns;
    fRefPatternMesh = cp.fRefPatternMesh;
    fNSubEl = cp.fNSubEl;
}

TPZRefPattern::TPZRefPattern(const TPZRefPattern &copy, const TPZPermutation &permute) :
        fRefPatternMesh(copy.fRefPatternMesh), fId(fNonInitializedId), fName("noname")
{
    fNSubEl = copy.fNSubEl;
    this->PermuteMesh(permute);

    SetRefPatternMeshToMasterDomain();

    ComputeTransforms();
    ComputePartition();
    GenerateSideRefPatterns();
}

TPZRefPattern::TPZRefPattern(TPZGeoMesh &gmesh){
    fRefPatternMesh = gmesh;
    CreateRefinementPattern();
}

TPZRefPattern::TPZRefPattern(std::istream &file) : fSideRefPattern(0), fId(fNonInitializedId), fName("noname")
{
    ReadAndCreateRefinementPattern(file);
}

TPZRefPattern::TPZRefPattern(const std::string &file ) : fSideRefPattern(0), fId(fNonInitializedId), fName("noname")
{
    std::ifstream input(file.c_str());
    ReadAndCreateRefinementPattern(input);

    /*
     BuildName();
     ofstream output(file.c_str());
     ExportPattern(output);
     */
}

int TPZRefPattern::operator==(const TPZAutoPointer<TPZRefPattern> compare) const
{
                if(fRefPatternMesh.NNodes() != compare->fRefPatternMesh.NNodes() || fRefPatternMesh.NElements() != compare->fRefPatternMesh.NElements())
                {
                                return 0;
                }
                TPZGeoEl *father = fRefPatternMesh.ElementVec()[0];
                TPZGeoEl *compfather = compare->fRefPatternMesh.ElementVec()[0];
                if(father->Type() != compfather->Type())
                {
                                return 0;
                }

                int nnodes = fRefPatternMesh.NNodes();
                int dim = father->Dimension();
                std::map<int,int> nodemap;
    //check coords
                REAL Tol;
                ZeroTolerance(Tol);
                int in;
                for(in = 0; in < nnodes; in++){
                                TPZManVector<REAL,3> coord(3,0.), coordcompare(3,0.);
                                TPZManVector<REAL,3> elparam(dim,0.), compareparam(dim,0.);

                                for(int i = 0; i < 3; i++){
                                                coord[i] = fRefPatternMesh.NodeVec()[in].Coord(i);
                                }

                                TPZManVector<REAL> diff(nnodes,0.);
                                int jn;
                                for(jn = 0; jn < nnodes; jn++){
                                                int j;
                                                for(j = 0; j < 3; j++){
                                                                coordcompare[j] = compare->fRefPatternMesh.NodeVec()[jn].Coord(j);
                                                }

                                                for(j=0 ; j<dim; j++){
                                                                diff[jn] += (coord[j]-coordcompare[j])*(coord[j]-coordcompare[j]);
                                                }

                                                diff[jn] = sqrt(diff[jn]);
                                                if(diff[jn] < ZeroTolerance()){
                                                                nodemap[in] = jn;
                                                                break;
                                                }
                                }

                                if(jn == nnodes)
                                {
                                                return 0;
                                }
                }

                int nelem = fRefPatternMesh.NElements();
                std::map<int,int> elementmap;

                int iel;
                for(iel = 0; iel < nelem; iel++)
                {
                                std::set<int> nodeset;
                                TPZGeoEl *igel = fRefPatternMesh.ElementVec()[iel];
                                int nnode = igel->NNodes();

                                int in;
                                for(in = 0; in < nnode; in++)
                                {
                                                nodeset.insert(nodemap[igel->NodeIndex(in)]);
                                }

                                int jel;
                                for(jel = 0; jel < nelem; jel++)
                                {
                                                if(elementmap.find(jel) != elementmap.end())
                                                {
                                                                continue;
                                                }
                                                std::set<int> compnodeset;
                                                TPZGeoEl *jgel = compare->fRefPatternMesh.ElementVec()[jel];
                                                int jnnode = jgel->NNodes();

                                                int jn;
                                                if(jnnode != nnode)
                                                {
                                                                continue;
                                                }

                                                for(jn = 0; jn < jnnode; jn++)
                                                {
                                                                compnodeset.insert(jgel->NodeIndex(jn));
                                                }

                                                if(nodeset == compnodeset)
                                                {
                                                                elementmap[jel] = iel;
                                                                break;
                                                }
                                }

                                if(jel == nelem)
                                {
                                                return 0;
                                }
                }

                return 1;
}


  void TPZRefPattern::Print(std::ostream &out) const
  {
     out << "TPZRefPattern::PrintMore\n\n";
                 int iSide,iSubEl;
                 out << "Refinement Pattern named " << fName << std::endl;

     out<<"SUB-ELEMENT INFO"<<std::endl;
     int nSubs = fSubElSideInfo.size();
                 for(iSubEl=0;iSubEl<nSubs;iSubEl++){
                     out<<"sub-el: "<<iSubEl<<std::endl;
                                 int nSides = fSubElSideInfo[iSubEl].size();
                                 for(iSide=0;iSide<nSides;iSide++){
             out<<"\tside: "<<iSide<<std::endl;
             out<<"\tfather side: "<<fSubElSideInfo[iSubEl][iSide].first<<std::endl;
                                     auto transform = fSubElSideInfo[iSubEl][iSide].second;
                                     out << "\tsub/side = " << iSubEl + 1 << "/" << iSide << "  FatherSide = " << fSubElSideInfo[iSubEl][iSide].first  << std::endl;
                                     out << "\tTransform = " << std::endl;
                                     transform.Mult().Print("Transformation T: ", out);
                                     transform.Sum().Print("Translation b: ", out);
                                 }
                 }

                 out<<"FATHER SIDES INFO"<<std::endl;
     int nSidesFather = fRefPatternMesh.ElementVec()[0]->NSides();
     for(iSide=0;iSide<nSidesFather;iSide++){
         out<<"side :"<<iSide<<std::endl;
         SPZFatherSideInfo &fatherInfo = fFatherSideInfo[iSide];

         out<<"\tInternal nodes:"<<std::endl;
         int nNodes = fatherInfo.fSideNodes.size();
         for(int iNode = 0; iNode < nNodes; iNode++){
             out<<"\t\tnode "<<iNode<<" index: "<<fatherInfo.fSideNodes[iNode]<<std::endl;
         }

         out<<"\tSubGeoElSides:"<<std::endl;
         int nSideSons = fatherInfo.fSideSons.size();
         for(int iGeoElSide = 0; iGeoElSide < nSideSons; iGeoElSide++){
            out<<"\t\telement: "<<FindSubEl(fatherInfo.fSideSons[iGeoElSide].Element()) + 1;
            out<<" side: "<< fatherInfo.fSideSons[iGeoElSide].Side()<<std::endl;
         }
     }
 }

void TPZRefPattern::PrintMore(std::ostream &out) const{
    Print(out);
    fRefPatternMesh.PrintTopologicalInfo(out);
}

void TPZRefPattern::ShortPrint(std::ostream &out) const
{
    TPZGeoEl *fatherEl = fRefPatternMesh.ElementVec()[0];
    int nsides = fatherEl->NSides();
    TPZVec<int> indices;
    //TPZVec<int> selected(nsides,0);
    out << fatherEl->TypeName();
    out << " Id " << fId << " Sides " ;
    for (int p=0 ; p<nsides; p++){
        if (fatherEl->SideDimension(p)==1 && NSideNodes(p))
        {
            out << p << " ";
        }
    }

    for(int n = 0; n < fRefPatternMesh.NodeVec().NElements(); n++)
    {
        out << " | ";
        for(int c = 0; c < 3; c++)
        {
            out << fRefPatternMesh.NodeVec()[n].Coord(c) << "   ";
        }
        out << " | ";
    }
}

void TPZRefPattern::PrintVTK(std::ofstream &file, bool matColor) const
{
    //          TPZGeoMesh * gmesh = &(RefPatternMesh());
    //          TPZRefPatternTools::PrintGMeshVTK(gmesh, file, matColor);

    TPZGeoMesh const & rGmesh = fRefPatternMesh;
    TPZGeoMesh * gmesh = new TPZGeoMesh;

    int64_t nNodes = rGmesh.NNodes();
    int64_t nElements = rGmesh.NElements();

    TPZVec < TPZVec <REAL> > NodeCoord(nNodes);
    for(int64_t i = 0; i < nNodes; i++) NodeCoord[i].Resize(3,0.);

    //setting nodes coords
    for(int64_t n = 0; n < nNodes; n++)
    {
        NodeCoord[n][0] = rGmesh.NodeVec()[n].Coord(0);
        NodeCoord[n][1] = rGmesh.NodeVec()[n].Coord(1);
        NodeCoord[n][2] = rGmesh.NodeVec()[n].Coord(2);
    }

    //initializing gmesh->NodeVec()
    gmesh->NodeVec().Resize(nNodes);
    TPZVec <TPZGeoNode> Node(nNodes);
    for(int64_t n = 0; n < nNodes; n++)
    {
        Node[n].SetNodeId(n);
        Node[n].SetCoord(NodeCoord[n]);
        gmesh->NodeVec()[n] = Node[n];
    }
    for(int64_t el = 1; el < nElements; el++)
    {
        TPZGeoEl * Elem = rGmesh.ElementVec()[el];
        TPZVec<int64_t> cornerindexes(Elem->NNodes());
        int64_t id = Elem->Id();
        for(int n = 0; n < Elem->NNodes(); n++)
        {
            cornerindexes[n] = Elem->NodeIndex(n);
        }
        gmesh->CreateGeoElement(Elem->Type(), cornerindexes, Elem->MaterialId(), id);
    }
    gmesh->BuildConnectivity();

    TPZVTKGeoMesh::PrintGMeshVTK(gmesh, file, matColor);
}

void TPZRefPattern::ExportPattern(std::ostream &out) const{
    out << fRefPatternMesh.NNodes() << ' '  << fRefPatternMesh.NElements() << std::endl << std::endl;
    out << fId << ' ' << fName << std::endl << std::endl ;

    for (int i = 0; i < fRefPatternMesh.NNodes() ; i++)
    {
        for (int k=0; k<3; k++)
        {
            out << fRefPatternMesh.NodeVec()[i].Coord(k) << ' ';
        }
        out << std::endl;
    }

    out << std::endl ;

    for (int i = 0; i < fRefPatternMesh.NElements(); i++)
    {
        out << fRefPatternMesh.ElementVec()[i]->Type() << " " << fRefPatternMesh.ElementVec()[i]->NCornerNodes() << " " ;
        for (int k = 0; k < fRefPatternMesh.ElementVec()[i]->NCornerNodes(); k++)
        {
            out << fRefPatternMesh.ElementVec()[i]->NodeIndex(k) << "  " ;
        }
        out << std::endl;
    }
    out << std::endl;
}

void TPZRefPattern::WritePattern(std::ofstream &out) const{

    ExportPattern(out);

    int nperm = fPermutedRefPatterns.size();
    out << nperm << ' ';
    for(int el = 0; el < nperm; el++)
    {
        out << fPermutedRefPatterns[el] << ' ';
    }
    out << std::endl;

    TPZGeoEl *father = fRefPatternMesh.ElementVec()[0];
    int nsides = father->NSides();
    int is;
    for(is=0; is<nsides; is++)
    {
        out << fSideRefPattern[is] << ' ';
    }
    out << std::endl;
}

int TPZRefPattern::ClassId() const {
    return Hash("TPZRefPattern");
}

void TPZRefPattern::Read(TPZStream &buf, void *context){
    int nSubElSides = -1;
    buf.Read(&nSubElSides);
    fSubElSideInfo.Resize(nSubElSides);
    for(int iSubElSide = 0; iSubElSide < nSubElSides; iSubElSide++){
        int nSides = -1;
        buf.Read(&nSides);
        fSubElSideInfo[iSubElSide].Resize(nSides);
        for(int iSide = 0; iSide < nSides; iSide++){
            int first = -1;
            buf.Read(&first);
            TPZTransform<REAL> second;
            second.Read(buf,context);
            fSubElSideInfo[iSubElSide][iSide] = std::make_pair(first,second);
        }
    }
    buf.Read(fFatherSideInfo);
    buf.Read(&fName);
    buf.Read(&fId);
    buf.Read(fSideRefPattern);
    buf.Read(fPermutedRefPatterns);
    fRefPatternMesh.Read(buf,context);
    buf.Read(&fNSubEl);
}

void TPZRefPattern::Write(TPZStream &buf, int withclassid) const{
    const int nSubElSides = fSubElSideInfo.size();
    buf.Write(&nSubElSides);
    for(int iSubElSide = 0; iSubElSide < nSubElSides; iSubElSide++){
        const int nSides = fSubElSideInfo[iSubElSide].size();
        buf.Write(&nSides);
        for(int iSide = 0; iSide < nSides; iSide++){
            buf.Write(&fSubElSideInfo[iSubElSide][iSide].first);
            fSubElSideInfo[iSubElSide][iSide].second.Write(buf,withclassid);
        }
    }
    buf.Write(fFatherSideInfo);
    buf.Write(&fName);
    buf.Write(&fId);
    buf.Write(fSideRefPattern);
    buf.Write(fPermutedRefPatterns);
    fRefPatternMesh.Write(buf,withclassid);
    buf.Write(&fNSubEl);

}


int TPZRefPattern::NSubElements() const
{
    #ifdef PZDEBUG
    if(fNSubEl != fSubElSideInfo.size()){
        PZError << "The refinement pattern is not in a consistent state."<<std::endl;
        PZError << "It has two different values of number of sub-elements, namely: "<<std::endl;
        PZError << fNSubEl<<" and "<<fSubElSideInfo.size()<<std::endl;
        PZError << "Aborting..."<<std::endl;
        DebugStop();
    }
    #endif
    return fNSubEl;
}

int TPZRefPattern::FatherSide(int side, int sub) const
{
    #ifdef PZDEBUG
    if(sub < 0 || sub >= NSubElements()){
        PZError << "TPZRefPattern::FatherSide: wrong parameter 'sub'."<<std::endl;
        PZError << "The current ref pattern has "<<NSubElements()<<" and the parameter 'sub' was "<<sub<<std::endl;
        DebugStop();

    }
    if(side < 0 || side >= fRefPatternMesh.ElementVec()[sub+1]->NSides()){
        PZError << "TPZRefPattern::FatherSide: wrong parameter 'side'."<<std::endl;
        PZError << "The current sub element has "<<fRefPatternMesh.ElementVec()[sub+1]->NSides()<<" sides ";
        PZError << "and the parameter 'side' was "<<side<<std::endl;
        DebugStop();
    }
    #endif
    return fSubElSideInfo[sub][side].first;
}
int TPZRefPattern::NSideSubGeoElSides(int fatherSide) const{
    #ifdef PZDEBUG
    if(!CheckSideConsistency(fatherSide)) {
        DebugStop();
    }
    #endif

    return fFatherSideInfo[fatherSide].fSideSons.size();
}

void TPZRefPattern::SideSubGeoElSide(int fatherSide, int subElPos, TPZGeoElSide & subGeoEl) const{
    #ifdef PZDEBUG
    if(!CheckSideAndSubElConsistency(fatherSide,subElPos)) {
        DebugStop();
    }
    #endif
    subGeoEl = fFatherSideInfo[fatherSide].fSideSons[subElPos];
}

 TPZTransform<> TPZRefPattern::Transform(int subElSide, int sub)
 {
    #ifdef PZDEBUG
     if(sub<0 || sub >= NSubElements()) {
         DebugStop();
     }
     TPZGeoEl *subEl = Element(sub+1);
     if(subElSide<0 || subElSide >= subEl->NSides()){
         DebugStop();
     }
    #endif
    return fSubElSideInfo[sub][subElSide].second;
 }

void TPZRefPattern::SideNodes(int fatherSide, TPZVec<int> &vecNodes){
#ifdef PZDEBUG
    if(!CheckSideConsistency(fatherSide)) {
        DebugStop();
    }
#endif
    vecNodes = fFatherSideInfo[fatherSide].fSideNodes;
}

int TPZRefPattern::NSideNodes(int fatherSide) const{
    #ifdef PZDEBUG
    if(!CheckSideConsistency(fatherSide)) {
        DebugStop();
    }
    #endif
    return fFatherSideInfo[fatherSide].fSideNodes.size();
}

int TPZRefPattern::NNodes() const{
    return fRefPatternMesh.NNodes();
}

bool TPZRefPattern::IsFatherNeighbour(TPZGeoElSide fatherSide,TPZGeoEl *son) const{
    const int side = fatherSide.Side();

    for (auto &geoElSide : fFatherSideInfo[side].fSideSons ) {
        if(geoElSide.Element() == son) return true;
    }
    return false;
}

TPZGeoEl *TPZRefPattern::Element(int iel){
    int nel = NSubElements()+1;/*sub-elements and father el*/
    if(iel < 0 || iel >= nel){
        PZError <<  "TPZRefPattern::Element the element with the following id does not exist: " << iel << std::endl;
        DebugStop();
    }
    return ( fRefPatternMesh.ElementVec()[iel]  );
}

void TPZRefPattern::InternalNodesIndexes(int side, TPZVec<TPZGeoElSideIndex> &nodeIndexes){
    #ifdef PZDEBUG
    if(!CheckSideConsistency(side)) {
        DebugStop();
    }
    #endif
    nodeIndexes.Resize(NSideNodes(side));
    int count = 0;
    for(auto &subElSide : fFatherSideInfo[side].fSideSons){
        if(subElSide.Dimension() == 0){
            nodeIndexes[count] = TPZGeoElSideIndex(subElSide);
            nodeIndexes[count].SetElementIndex(nodeIndexes[count].ElementIndex());
            count++;
        }
    }
}

void TPZRefPattern::InternalSidesIndexes(int side, TPZVec<TPZGeoElSideIndex> &sideIndexes){
    #ifdef PZDEBUG
    if(!CheckSideConsistency(side)) {
        DebugStop();
    }
    #endif
    const int nSubElSides = fFatherSideInfo[side].fSideSons.size();
    sideIndexes.Resize(nSubElSides);
    for(int iSub = 0; iSub < nSubElSides; iSub++){
        sideIndexes[iSub] = TPZGeoElSideIndex(fFatherSideInfo[side].fSideSons[iSub]);
        sideIndexes[iSub].SetElementIndex(sideIndexes[iSub].ElementIndex());
    }
}

int TPZRefPattern::SidePartition(TPZVec<TPZGeoElSide> &gelvec, int side)
{
    #ifdef PZDEBUG
    if(!CheckSideConsistency(side)) {
        DebugStop();
    }
    #endif
    gelvec = fFatherSideInfo[side].fSideSons;
    return gelvec.size();//numero de elementos da particao do lado
}

TPZAutoPointer<TPZRefPattern> TPZRefPattern::SideRefPattern(int side)
{
                const int id = this->fSideRefPattern[side];

                return gRefDBase.FindRefPattern(id);
}

void TPZRefPattern::CreateNewNodes(TPZGeoEl * gel, TPZVec<int64_t> &newnodeindexes)
{
                int nNodes = gel->NCornerNodes();
                int totalNodes = fRefPatternMesh.NNodes();
                newnodeindexes.Resize(totalNodes);

                if (gel->HasSubElement()){
                                std::cout << "CreateNewNodes called for an element which is already divided. \n";
                                std::cout.flush();
                                return;
                }

                int nSides = gel->NSides();
                for (int side = nNodes;side<nSides;side++){
                                CreateMidSideNodes(gel,side,newnodeindexes);
                }
}

void TPZRefPattern::CreateMidSideNodes(TPZGeoEl * gel, int side, TPZVec<int64_t> &newnodeindexes)
{
                TPZGeoMesh *gmesh = gel->Mesh();
                //SideNodes returns a vector with the indexes of the internal nodes of the side side
                TPZManVector<int> sideNodes;
                SideNodes(side,sideNodes);
                TPZGeoElSide gelside(gel,side);
                TPZGeoElSide neighbour(gelside.Neighbour());
                TPZManVector<int64_t> sideindices(0);
                //checks if a neighbour element has created a midsidenode already
                while(neighbour.Element() && neighbour != gelside){
                                if(neighbour.HasSubElement() && neighbour.Element()->NSideSubElements(neighbour.Side()) > 1)
                                {
                                                neighbour.Element()->MidSideNodeIndices(neighbour.Side(),sideindices);
                                                break;
                                }
                                neighbour = neighbour.Neighbour();
                }
                for (int j=0;j<sideNodes.NElements();j++){
                                int index = sideNodes[j];
                                //new node coordinates
                                TPZManVector<REAL,3> refnodecoord_X(3,0.);
                                TPZManVector<REAL,3> neighbourcoord_X(3,0.);

                                for (int k=0;k<3;k++){
            double coord = fRefPatternMesh.NodeVec()[index].Coord(k);
            refnodecoord_X[k] = coord;
        }

                                TPZManVector<REAL,3> newnodecoord_xi(Element(0)->Dimension(),0.);

                                //getting node coordinates in refpattern mesh: just a matter of adjusting dimension
                                for(int iX = 0; iX < Element(0)->Dimension(); iX++) newnodecoord_xi[iX] = refnodecoord_X[iX];

        //coordinates of the node to be created (if it hasnt been created yet)
                                TPZManVector<REAL,3> newnodecoord_X(3);
                                gel->X(newnodecoord_xi,newnodecoord_X);

                                newnodeindexes[index] = -1;
                                REAL smallestDiff = -1.;
                                int smallestDiffIndex = -1;
                                //check if the node has already been created using information previsouly obtained from neighbour
                                for(int i=0; i< sideindices.NElements(); i++)
                                {
                                                for(int k=0; k<3; k++) neighbourcoord_X[k] = gmesh->NodeVec()[sideindices[i]].Coord(k);
                                                REAL dif = 0.;
                                                for (int k=0;k<3;k++)
                                                {
                                                                dif += (newnodecoord_X[k] - neighbourcoord_X[k]) * (newnodecoord_X[k] - neighbourcoord_X[k]);
                                                }
                                                if(smallestDiffIndex < 0. || smallestDiff > dif)
                                                {
                                                                smallestDiff = dif;
                                                                smallestDiffIndex = i;
                                                }
                                }
                                if (smallestDiff < 1e-2 && sideindices.NElements() != 0)
                                {
                                                newnodeindexes[index] = sideindices[smallestDiffIndex];
                                }
                                if (smallestDiff >= 1.e-2 && sideindices.NElements() != 0)
                                {
#ifdef LOG4CXX
                                                {
                                                                std::stringstream sout;
                                                                sout << "Incompatible refinement patterns detected\n";
                                                                sout << "Closest node at distance " << smallestDiff << std::endl;
                                                                gel->Print(sout);
                                                                LOGPZ_ERROR(logger,sout.str())
                                                }
#endif
            std::cout << "Refpattern is trying to create midnode but there is another node which is not in the right position!\n";
            std::cout << "Either the refinement patterns are incompatible or there is an issue with the mapping.\n";
            std::cout << "newnodecoord_X(" << newnodecoord_X[0] << "," << newnodecoord_X[1] << "," << newnodecoord_X[2] << ")" << std::endl;
            std::cout << "neighbourcoord_X(" << neighbourcoord_X[0] << "," << neighbourcoord_X[1] << "," << neighbourcoord_X[2] << ")" << std::endl;
                                                DebugStop();
                                }
                                if (newnodeindexes[index] == -1)
                                {
                                                //If no one has created the node yet, it should be created
                                                int64_t newindex = gmesh->NodeVec().AllocateNewElement();
                                                gmesh->NodeVec()[newindex].Initialize(newnodecoord_X,*gmesh);
                                                newnodeindexes[index] = newindex;
                                }
                }
}


void TPZRefPattern::GenerateSideRefPatterns(){
    if(!fRefPatternMesh.ElementVec().NElements())
    {
        return;
    }

    TPZGeoEl *fatherEl = fRefPatternMesh.ElementVec()[0];

    int nsides = fatherEl->NSides();
    fSideRefPattern.Resize(nsides,-1);
    fSideRefPattern.Fill(-1.);

    int thisId = this->Id();
    fSideRefPattern[nsides-1] = thisId;

    int is;
    for(is=0; is<nsides-1; is++)
    {
        if(fatherEl->SideDimension(is) == 0)  continue;
        if(NSideSubGeoElSides(is) == 1) continue;

        TPZAutoPointer<TPZRefPattern> sideref(new TPZRefPattern());

        BuildSideMesh(is, sideref->fRefPatternMesh);

        sideref->SetRefPatternMeshToMasterDomain();

        sideref->fNSubEl = sideref->fRefPatternMesh.NElements()-1;
        sideref->ComputeTransforms();
        sideref->ComputePartition();
        TPZAutoPointer<TPZRefPattern> found = gRefDBase.FindRefPattern(sideref);

        if(!found.operator ->())
        {
            if(this->NameInitialized())
            {
                sideref->BuildName();
            }
            sideref->GenerateSideRefPatterns();
            gRefDBase.InsertRefPattern(sideref);
            fSideRefPattern[is] = sideref->Id();
            sideref->InsertPermuted();
        }
        else
        {
            fSideRefPattern[is] = found->Id();
        }
    }
}

void TPZRefPattern::InsertPermuted()
{
    if(!fRefPatternMesh.ElementVec().NElements() || !fRefPatternMesh.ElementVec()[0])
    {
        return;
    }

    GenerateSideRefPatterns();
    TPZGeoEl * gel = fRefPatternMesh.ElementVec()[0];
    GeneratePermutations(gel);
    MElementType geltype = gel->Type();

    std::list<TPZRefPatternPermute> &permlist = fPermutations[geltype];
    std::list<TPZRefPatternPermute>::iterator it;
    fPermutedRefPatterns.resize(permlist.size());

    int64_t counter;
    for(it=permlist.begin(), counter=0; it != permlist.end(); it++,counter++)
    {
        TPZAutoPointer<TPZRefPattern> refp(new TPZRefPattern(*this,(*it).fPermute));
        TPZAutoPointer<TPZRefPattern> found = gRefDBase.FindRefPattern(refp);

#ifdef LOG4CXX
        if (logger->isDebugEnabled())
        {
            std::stringstream sout;
            sout << "Permutation " << it->fPermute;
            sout << "Created refpattenr refp ";
            refp->Print(sout);
            sout << "found refpattern ";
            if(found) found->Print(sout);
            else sout << "Pattern not found";
            LOGPZ_DEBUG(logger,sout.str())
        }
#endif

        if(found)
        {
            fPermutedRefPatterns[counter] = found->Id();
        }
        else
        {
            if(this->NameInitialized())
            {
                refp->BuildName();
            }
            gRefDBase.InsertRefPattern(refp);
        }
    }
}

TPZAutoPointer<TPZRefPattern> TPZRefPattern::SideRefPattern(int side, TPZTransform<> &trans){
    TPZAutoPointer<TPZRefPattern> sideref = this->SideRefPattern(side);
    if(!sideref) return nullptr;
    return sideref->FindRefPattern(trans);
}

/****************************PROTECTED METHODS*************************************************************************/

TPZAutoPointer<TPZRefPattern> TPZRefPattern::FindRefPattern(TPZTransform<> &trans)
{
    REAL tol = 1.e-6;
    if(!fRefPatternMesh.ElementVec().NElements() || ! fRefPatternMesh.ElementVec()[0]) return 0;
    MElementType type = fRefPatternMesh.ElementVec()[0]->Type();

    std::list<TPZRefPatternPermute> &permlist = fPermutations[type];
    std::list<TPZRefPatternPermute>::iterator it;
    for(it = permlist.begin(); it != permlist.end(); it++)
    {
        TPZRefPatternPermute &tmp = (*it);
        if(!tmp.fTransform.CompareTransform(trans, tol))
        {
            TPZAutoPointer<TPZRefPattern> refpPerm = new TPZRefPattern(*this, tmp.fPermute);
            TPZAutoPointer<TPZRefPattern> found = gRefDBase.FindRefPattern(refpPerm);

            return found;
        }
    }
    return 0;
}

void TPZRefPattern::SetRefPatternMeshToMasterDomain(){
    TPZGeoEl * fatherEl = fRefPatternMesh.ElementVec()[0];


    const int dim = fatherEl->Dimension();
    const int nNodes = fRefPatternMesh.NodeVec().NElements();
    TPZManVector< TPZManVector<REAL,3> ,20> nodecoords_inQSI(nNodes,TPZManVector<REAL,3>(dim,0));
    TPZManVector<REAL,3> nodecoords_inX(3,0);
    TPZManVector<REAL,3> temp(3,0.);
    REAL Tol;
    ZeroTolerance(Tol);
    for(int n = 0; n < nNodes; n++)
    {
        fRefPatternMesh.NodeVec()[n].GetCoordinates(nodecoords_inX);
        fatherEl->ComputeXInverse(nodecoords_inX, nodecoords_inQSI[n],Tol);
    }

    TPZVec<REAL> coordQSIprojected(dim,0.);
    for(int n = 0; n < nNodes; n++)
    {
        //Let us guarantee that ComputeXInverse did not result in a point outside the reference element
        fatherEl->ProjectInParametricDomain(nodecoords_inQSI[n], coordQSIprojected);

        int c;
        for(c = 0; c < dim; c++)
        {
            double qsi = coordQSIprojected[c];
            fRefPatternMesh.NodeVec()[n].SetCoord(c,qsi);
        }
        for(; c < 3; c++)
        {
            fRefPatternMesh.NodeVec()[n].SetCoord(c,0.);
        }
    }
}

void TPZRefPattern::GeneratePermutations(TPZGeoEl *fatherEl){
    if (!fatherEl)
    {
        PZError << "Error at " << __PRETTY_FUNCTION__ << " at line " << __LINE__ << " Father Element is NULL\n";
        return;
    }
    MElementType fatherEltype = fatherEl->Type();

    if(fPermutations.count(fatherEltype))
    {
        return;
    }
    TPZGeoMesh *fatherElMesh = fatherEl->Mesh();
    TPZGeoMesh gmesh;
    gmesh.NodeVec().Resize(fatherEl->NNodes());

    int iNode, nNodes = fatherEl->NNodes();
    for(iNode = 0; iNode < nNodes; iNode++)
    {
        gmesh.NodeVec()[iNode].Initialize(fatherElMesh->NodeVec()[fatherEl->NodeIndex(iNode)], gmesh);
    }
    TPZManVector<int64_t,8> nodes(nNodes), nodesPerm(nNodes);
    for(iNode = 0; iNode < nNodes; iNode++)
    {
        nodes[iNode]=iNode;
    }
    int matid = fatherEl->MaterialId();
    TPZPermutation permute(nNodes);
    TPZTransform<> trans(fatherEl->Dimension());
    TPZRefPatternPermute refpermute;

    refpermute.fPermute = permute;
    refpermute.fTransform = trans;
    fPermutations[fatherEl->Type()].push_back(refpermute);
    permute++;
    TPZGeoEl *fatherElPermuted;
    TPZVec<int> permutedSides;
    while(!permute.IsFirst())
    {
        permute.Permute(nodes,nodesPerm);
        int64_t index;
        fatherElPermuted = gmesh.CreateGeoElement(fatherEl->Type(),nodesPerm,matid,index,0);
//        fatherEl->GetPermutation()
        bool valid = false;
        int nPermutes = fatherEl->NPermutations();
        for(int iPermute = 1; iPermute < nPermutes; iPermute++){//its ok to skip the first (trivial) permutation
            fatherEl->GetPermutation(iPermute,permutedSides);
            for(iNode = 0; iNode < nNodes; iNode++) {
                if ( permutedSides[iNode] != nodesPerm[iNode]) break;
            }
            if(iNode == nNodes){
                valid = true;//found the correspondent permutation
                break;
            }
        }
        if(valid)
        {
            TPZRefPatternPermute candidate;
            candidate.fPermute = permute;
            candidate.fTransform = ComputeParamTransform(fatherEl,fatherElPermuted,fatherEl->NSides()-1,fatherEl->NSides()-1);
            fPermutations[fatherEl->Type()].push_back(candidate);
        }
        permute++;
    }
}

void TPZRefPattern::ComputeTransforms(){
    TPZGeoEl *fatherEl = fRefPatternMesh.Element(0);
    if(!fatherEl)
    {
        PZError << "TPZRefPattern::ComputePartition Father not exists?!\n";
        DebugStop();
    }

    const int nSubEls = fRefPatternMesh.NElements() - 1;
    fSubElSideInfo.Resize(nSubEls);

    REAL Tol;
    ZeroTolerance(Tol);
    const int dim = fatherEl->Dimension();
    TPZManVector<REAL,3> massCenter(3,0.), xPoint(3,0.), fatherParam(dim,0.);
    for(int iSubEl=0; iSubEl < nSubEls; iSubEl++)
    {
        TPZGeoEl *son = Element(iSubEl+1);
        massCenter.Resize(son->Dimension());
        int nSides = son->NSides();
        fSubElSideInfo[iSubEl].Resize(nSides);
        for(int iSide = 0; iSide < nSides; iSide++)
        {
            TPZGeoElSide elside (son, iSide);

            son->CenterPoint(iSide,massCenter);
            son->X(massCenter,xPoint);//already in father's coordinates. but the dimension must be set up
            for(int iX = 0; iX < dim; iX++) fatherParam[iX] = xPoint[iX];
            const int fatherSide = fatherEl->WhichSide(fatherParam);
            fSubElSideInfo[iSubEl][iSide].first = fatherSide;
//            fSubElSideInfo[iSubEl][iSide].second = son->ComputeParamTrans(fatherEl,fatherSide,iSide);
            fSubElSideInfo[iSubEl][iSide].second = ComputeParamTransform(son,fatherEl,iSide,fatherSide);
        }
    }
}

void TPZRefPattern::ComputePartition(){
    //this method will fill the fFatherSideInfo data structure
    TPZGeoEl *fatherEl = Element(0);
    int nSides = fatherEl->NSides();

    fFatherSideInfo.Resize(nSides);
    for(int iSide = 0; iSide < nSides; iSide++){

        TPZStack<TPZGeoElSide> sideSonsStack;
        TPZStack<int> cornerNodeStack;
        TPZStack<int> internalNodesStack;
        const int nSideNodes = fatherEl->NSideNodes(iSide);
        for(int iNode = 0; iNode < nSideNodes; iNode++){
            const int nodeIndex = fatherEl->SideNodeIndex(iSide,iNode);
            cornerNodeStack.push_back(nodeIndex);
        }
        for(auto iSubEl = 0; iSubEl < NSubElements(); iSubEl++){
            TPZGeoEl *subEl = Element(iSubEl+1);
            for(auto iSubElSide = 0; iSubElSide < subEl->NSides(); iSubElSide ++){
                int fatherSide = FatherSide(iSubElSide,iSubEl);
                if(fatherSide == iSide){
                    TPZGeoElSide geoElSideCandidate(subEl,iSubElSide);
                    bool isInList = false;
                    for(int iGelSide = 0; iGelSide < sideSonsStack.size(); iGelSide++){
                        if(isInList) break;
                        TPZGeoElSide neighbour = geoElSideCandidate.Neighbour();
                        while(neighbour.Element() && neighbour.Element()->Id() != geoElSideCandidate.Element()->Id()){
                            if(neighbour == sideSonsStack[iGelSide]){
                                isInList = true;
                                break;
                            }
                            neighbour = neighbour.Neighbour();
                        }
                    }
                    if(!isInList)   {
                        sideSonsStack.push_back(geoElSideCandidate);
                        if(subEl->SideDimension(iSubElSide) == 0){
                            const int nodeIndex = subEl->NodeIndex(iSubElSide);
                            bool isNodeInList = false;
                            const int nCornerNodes = cornerNodeStack.size();
                            const int nInternalNodes = internalNodesStack.size();
                            for(int iNode = 0; iNode < nInternalNodes; iNode++){
                                if(internalNodesStack[iNode] == nodeIndex) isNodeInList = true;
                            }
                            //the following garantees that for a 0D side its node is put correctly in internalNodesStack
                            for(int iNode = 0; iNode < nCornerNodes && fatherEl->SideDimension(iSide) > 0; iNode++){
                                if(cornerNodeStack[iNode] == nodeIndex) isNodeInList = true;
                            }
                            if(!isNodeInList){
                                internalNodesStack.push_back(nodeIndex);
                            }
                        }
                    }
                }
            }
        }
        fFatherSideInfo[iSide].fSideNodes.Resize(internalNodesStack.size());
        for(auto iNode = 0; iNode < internalNodesStack.size(); iNode++){
            fFatherSideInfo[iSide].fSideNodes[iNode] = internalNodesStack[iNode];
        }
        fFatherSideInfo[iSide].fSideSons.Resize(sideSonsStack.size());
        for(auto iSubEl = 0; iSubEl < sideSonsStack.size(); iSubEl++){
            fFatherSideInfo[iSide].fSideSons[iSubEl] = sideSonsStack[iSubEl];
        }
    }
}

 int TPZRefPattern::FindSubEl(TPZGeoEl *geoEl) const{
    const int nSubEls = NSubElements();
    for(int i = 0; i < nSubEls; i++){
        if(fRefPatternMesh.ElementVec()[i+1] == geoEl){
            return i;
        }
    }
    PZError<<"TPZRefPattern::FindSubEl: ERROR"<<std::endl;
    PZError<<"This method should not be called with an element ";
    PZError<<"that is not a sub-element of this Refinement Pattern."<<std::endl;
    PZError<<"Aborting..."<<std::endl;
    DebugStop();
    return -1;
}

bool TPZRefPattern::CheckSideConsistency(const int fatherSide) const{
    int nSides = fRefPatternMesh.ElementVec()[0]->NSides();
    if(fatherSide < 0 || fatherSide >= nSides)
    {
        PZError << "TPZRefPattern::NSideSubElements: wrong argument\n";
        PZError << "father side = " << fatherSide << std::endl;
        return false;
    }
    return true;
}

bool TPZRefPattern::CheckSideAndSubElConsistency(const int fatherSide, const int subEl) const{
    const int nSides = fRefPatternMesh.ElementVec()[0]->NSides();
    const int nSubGeoElSides = NSideSubGeoElSides(fatherSide);
    if(fatherSide < 0 || fatherSide >= nSides || subEl < 0 || subEl >= nSubGeoElSides)
    {
        PZError << "TPZRefPattern::CheckSideAndSubElConsistency: wrong argument\n";
        PZError << "father side = " << fatherSide << std::endl;
        PZError << "subEl = " << subEl << std::endl;
        return false;
    }
    return true;
}

void TPZRefPattern::BuildSideMesh(int side, TPZGeoMesh &SideRefPatternMesh)
{
    if(!fRefPatternMesh.ElementVec().NElements())
    {
        return;
    }
    TPZGeoEl *gel = fRefPatternMesh.ElementVec()[0];
    TPZStack<int> allsides;
    std::map<int,int> allsidenodes;
    int count =0;
    gel->LowerDimensionSides(side, allsides);
    allsides.Push(side);
    int s;
    for(s=0; s<allsides.NElements(); s++)
    {
        TPZStack<int> sidenodes;
        SideNodes(allsides[s],sidenodes);
        int t;
        for(t=0; t<sidenodes.NElements(); t++)
        {
            if(!allsidenodes.count(sidenodes[t])) allsidenodes[sidenodes[t]] = count++;
        }
    }
    SideRefPatternMesh.NodeVec().Resize(allsidenodes.size());
    std::map<int,int>::iterator it;
    for(it = allsidenodes.begin(); it!= allsidenodes.end(); it++)
    {
        int nodeorig = (*it).first;
        int nodedest = (*it).second;

        SideRefPatternMesh.NodeVec()[nodedest].Initialize(fRefPatternMesh.NodeVec()[nodeorig], SideRefPatternMesh);
        SideRefPatternMesh.NodeVec()[nodedest].SetNodeId(nodedest);
    }
    TPZStack<int64_t> nodeindices;
    nodeindices.Resize(gel->NSideNodes(side));
    int64_t in;
    for(in=0; in<nodeindices.NElements(); in++)
    {
        nodeindices[in] = allsidenodes[gel->SideNodeIndex(side,in)];
    }
    int matid = gel->MaterialId();
    int64_t index;
    TPZGeoEl *father = SideRefPatternMesh.CreateGeoElement(gel->Type(side),nodeindices,matid,index,1);
    int sidedim = father->Dimension();
    TPZStack<TPZGeoElSide> gelvec;
    SidePartition(gelvec, side);
    int nsub = gelvec.NElements();
    int subel;
    for(subel=0; subel<nsub; subel++)
    {
        if(gelvec[subel].Dimension() != sidedim)
        {
            continue;
        }
        nodeindices.Resize(gelvec[subel].NSideNodes());
        for(in=0; in<nodeindices.NElements(); in++)
        {
            nodeindices[in] = allsidenodes[gelvec[subel].SideNodeIndex(in)];
        }
        MElementType type = gelvec[subel].Element()->Type(gelvec[subel].Side());
        TPZGeoEl *subel = SideRefPatternMesh.CreateGeoElement(type,nodeindices,matid,index);
        subel->SetFather(father);
    }
    SideRefPatternMesh.BuildConnectivity();
}


void TPZRefPattern::BuildName()
{
    fName = TPZRefPatternTools::BuildRefPatternModelName(*this);
}

void TPZRefPattern::PermuteMesh(const TPZPermutation &permute)
{
    int in;
    TPZGeoEl *father = this->fRefPatternMesh.ElementVec()[0];
    int nn = father->NNodes();
    TPZVec<int> nodes(nn), nodeperm(nn);
    for(in=0; in<nn; in++){
        nodes[in] = father->NodeIndex(in);
    }
    permute.Permute(nodes,nodeperm);
    for(in=0; in<nn; in++){
        father->SetNodeIndex(in, nodeperm[in]);
    }
    fRefPatternMesh.ResetConnectivities();
    fRefPatternMesh.BuildConnectivity();
}

void TPZRefPattern::ReadAndCreateRefinementPattern(std::istream &pattern){
    int nNodes, nElems;
    pattern >> nNodes >> nElems;
    pattern >> fId >> fName;

    TPZManVector<REAL,3> coord(3);
    fRefPatternMesh.NodeVec().Resize(nNodes);

    //criacao dos nohs
    for(int iNode = 0; iNode < nNodes; iNode++)
    {
        pattern >> coord[0];
        pattern >> coord[1];
        pattern >> coord[2];
        fRefPatternMesh.NodeVec()[iNode].Initialize(iNode,coord,fRefPatternMesh);
    }
    TPZGeoEl *father = 0;
    //criacao dos elementos geometricos que definem a particao
    int ntype, nummat, naorners, incid, el;
    for(el=0; el<nElems; el++)//os sub-elementos podem nao ter de uma mesma geometria
    {
        pattern >> ntype >> nummat;
        MElementType etype = (MElementType) ntype;
        naorners = MElementType_NNodes(etype);
        TPZVec<int64_t> nodes(naorners);
        for(incid = 0; incid < naorners; incid++)
        {
            pattern >> nodes[incid];
        }
        int64_t index;
        TPZGeoEl *subel = fRefPatternMesh.CreateGeoElement(etype, nodes, nummat, index, 0);
        if(el == 0)
        {
            father = subel;
            fRefPatternMesh.SetDimension(father->Dimension());
        }
        if(el > 0)
        {
            subel->SetFather(father);
            subel->SetFatherIndex(father->Index());
        }
    }
    CreateRefinementPattern();
}

void TPZRefPattern::CreateRefinementPattern(){
    //ensure that the all the sub-elements are aware of their father
    //the first element is supposed to be the father element
    TPZGeoEl* fatherEl = fRefPatternMesh.Element(0);
    MElementType type = fatherEl->Type();
    fNSubEl = fRefPatternMesh.NElements() - 1;

    for(int iSubEl = 0; iSubEl < fNSubEl; iSubEl++){
        TPZGeoEl* subEl = fRefPatternMesh.Element(iSubEl + 1);
        subEl->SetFather(fatherEl);
        subEl->SetFatherIndex(0);
    }

    SetRefPatternMeshToMasterDomain();
    fRefPatternMesh.BuildConnectivity();//sub-element connectivities
    GeneratePermutations(fatherEl);

    ComputeTransforms();//compute the transformation between sons and father-el
    ComputePartition();//compute the partition of the father element in (side,sub-element)

    GenerateSideRefPatterns();
}

TPZTransform<> TPZRefPattern::ComputeParamTransform(TPZGeoEl *geoElSon, TPZGeoEl *geoElFather, int sideSon, int sideFather){
    int dimFather = geoElFather->Dimension();
    int dimSideFather = geoElFather->SideDimension(sideFather);
    int dimSon = geoElSon->Dimension();
    int dimSideSon = geoElSon->SideDimension(sideSon);
    if(dimSideFather < dimSideSon){
        PZError << "\nTPZRefPattern::ComputeParamTransform called with sides error\n";
        DebugStop();
    }

    if(!geoElFather->SideDimension(sideFather)) return TPZTransform<>(0,0);

    REAL weight;
    TPZFNMatrix<9> jac(dimSon,dimSon),axes(3,3,0.);
    TPZFNMatrix<9> jacinv(dimSon,dimSon);
    TPZManVector<REAL,3> x(3,0.);
    TPZManVector<REAL,3> intpoint(dimSideSon,0.);
    int tam = (dimSideSon+1);
    TPZFNMatrix<16> hess(tam,tam,0.),grad0(tam,1,0.);
    TPZIntPoints *intrule = geoElSon->CreateSideIntegrationRule(sideSon,2);
    TPZManVector<int,3> order(dimSideSon,2);
    intrule->SetOrder(order);
    //integrating over the side of the son that is contained on the father's side
    int ij,ik,indp;
    REAL D2Edaikdaij,D2Edcidaij,D2Edci2;
    D2Edcidaij = 0.;
    D2Edci2 = 0.;
    for(ij=0;ij<dimSideSon;ij++){
        for(ik=ij;ik<dimSideSon;ik++){
            D2Edaikdaij = 0.;
            D2Edcidaij  = 0.;
            for(indp = 0; indp < intrule->NPoints(); ++indp){
                intrule->Point(indp,intpoint,weight);
                D2Edaikdaij += intpoint[ik]*intpoint[ij]*weight;
                if(ik==ij) D2Edcidaij += intpoint[ij]*weight;
                if(ij==0 && ik==0) D2Edci2 += weight;
            }
            hess(ij,ik) = 2.*D2Edaikdaij;
            hess(ik,ij) = hess(ij,ik);
            if(ik==ij) {
                hess(ij,dimSideSon) = 2.*D2Edcidaij;
                hess(dimSideSon,ij) = hess(ij,dimSideSon);
            }
            if(ij==0 && ik==0) hess(dimSideSon,dimSideSon) = 2.*D2Edci2;
        }
    }// end of integral
    //defining transform from the son's side to the son's element
    TPZTransform<> transformSideToSon(dimSon);//identidade
    if(dimSideSon<dimSon) transformSideToSon = geoElSon->SideToSideTransform(sideSon,geoElSon->NSides()-1);
    TPZTransform<> fatelside = geoElFather->SideToSideTransform(geoElFather->NSides()-1,sideFather);
    TPZManVector<REAL,3> sidepoint(dimSon);
    int j;
    TPZFNMatrix<9> Amat(dimSideFather,dimSideSon,0.),Bmat(dimSideFather,1,0.);
    REAL Tol;
    ZeroTolerance(Tol);
    TPZVec<REAL> xiFather(dimFather,0.);
    TPZVec<REAL> xiFatherSide(dimSideFather,0.);
    for(int ifat=0;ifat<dimSideFather;ifat++){
        REAL DEdci = 0.;
        for(j=0;j<(dimSideSon+1);j++){
            REAL DEdaij = 0.;
            for(indp = 0; indp < intrule->NPoints(); ++indp){
                intrule->Point(indp,intpoint,weight);
                transformSideToSon.Apply(intpoint,sidepoint);//lado do filho para o seu interior: mestre
                geoElSon->X(sidepoint,x);//this is already on father's domain, but in 3 dimensions
                for(int iX = 0; iX < dimFather; iX++) xiFather[iX] = x[iX];
                fatelside.Apply(xiFather,xiFatherSide);
                if(j<dimSideSon) DEdaij += xiFatherSide[ifat]*intpoint[j]*weight;
                if(j==0) DEdci += xiFatherSide[ifat]*weight;
            }
            if(j<dimSideSon) grad0(j,0) = 2.*DEdaij;
            if(j==0) grad0(dimSideSon,0) = 2.*DEdci;
            if(!dimSideSon) grad0(0,0) = DEdci;
        }

        if(dimSideSon) hess.SolveDirect(grad0,ELU);
        for(int k=0;k<dimSideSon;k++) Amat(ifat,k) = grad0(k,0);
        Bmat(ifat,0) = grad0(dimSideSon,0);
    }

    delete intrule;
    TPZTransform<> t(dimSideFather,dimSideSon);
    t.SetMatrix(Amat,Bmat);
    return t;
}

int TPZRefPattern::TPZRefPatternPermute::ClassId() const {
 return Hash("TPZRefPattern::TPZRefPatternPermute");
}

void TPZRefPattern::TPZRefPatternPermute::Read(TPZStream& buf, void* context) { //ok
 fPermute.Read(buf, context);
 fTransform.Read(buf, context);
}

void TPZRefPattern::TPZRefPatternPermute::Write(TPZStream& buf, int withclassid) const { //ok
 fPermute.Write(buf, withclassid);
 fTransform.Write(buf, withclassid);
}

void TPZRefPattern::SPZFatherSideInfo::Read(TPZStream& buf, void* context) { //ok
    buf.Read(fSideNodes);
    buf.Read(fSideSons);
}

void TPZRefPattern::SPZFatherSideInfo::Write(TPZStream& buf, int withclassid) const { //ok
    buf.Write(fSideNodes);
    buf.Write(fSideSons);
}