/[escript]/trunk/finley/src/CPPAdapter/MeshAdapter.cpp
ViewVC logotype

Diff of /trunk/finley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1361 by gross, Fri Dec 14 09:26:51 2007 UTC revision 1375 by gross, Wed Jan 9 00:15:05 2008 UTC
# Line 809  void MeshAdapter::addPDEToRHS( escript:: Line 809  void MeshAdapter::addPDEToRHS( escript::
809     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );
810     checkFinleyError();     checkFinleyError();
811  }  }
812    //
813    // adds PDE of second order into a transport problem
814    //
815    void MeshAdapter::addPDEToTransportProblem(
816                         TransportProblemAdapter& tp, escript::Data& source, const escript::Data& M,
817                         const escript::Data& A, const escript::Data& B, const escript::Data& C,const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
818                         const escript::Data& d, const escript::Data& y,
819                         const escript::Data& d_contact,const escript::Data& y_contact) const
820    {
821       DataArrayView::ShapeType shape;
822       source.expand();
823       escript:: Data tmp(0.0,M.getDataPointShape(),tp.getFunctionSpace(),true);
824       escriptDataC _source=source.getDataC();
825       escriptDataC _tmp=tmp.getDataC();
826       escriptDataC _M=M.getDataC();
827       escriptDataC _A=A.getDataC();
828       escriptDataC _B=B.getDataC();
829       escriptDataC _C=C.getDataC();
830       escriptDataC _D=D.getDataC();
831       escriptDataC _X=X.getDataC();
832       escriptDataC _Y=Y.getDataC();
833       escriptDataC _d=d.getDataC();
834       escriptDataC _y=y.getDataC();
835       escriptDataC _d_contact=d_contact.getDataC();
836       escriptDataC _y_contact=y_contact.getDataC();
837    
838       Finley_Mesh* mesh=m_finleyMesh.get();
839       Paso_FCTransportProblem* _tp = tp.getPaso_FCTransportProblem();
840      
841    
842       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_tmp, &_M);
843       checkFinleyError();
844       /* add mass matix to lumped mass matrix of transport problem */
845       double* tmp_prt=getSampleData(&_tmp,0);
846       int i;
847       int n=Paso_FCTransportProblem_getTotalNumRows(_tp);
848       #pragma omp parallel for private(i) schedule(static)
849       for (i=0;i<n ;++i) _tp->lumped_mass_matrix[i]+=tmp_prt[i];
850    
851       Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->transport_matrix, &_source, &_A, 0, 0, &_D, &_X, &_Y );
852       checkFinleyError();
853    
854       Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->flux_matrix, &_source, 0, &_B, &_C, 0, 0, 0 );
855       checkFinleyError();
856    
857       Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, _tp->transport_matrix, &_source, 0, 0, 0, &_d, 0, &_y );
858       checkFinleyError();
859    
860       Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );
861       checkFinleyError();
862    }
863    
864  //  //
865  // interpolates data between different function spaces:  // interpolates data between different function spaces:
# Line 920  void MeshAdapter::interpolateOnDomain(es Line 971  void MeshAdapter::interpolateOnDomain(es
971          }          }
972          break;          break;
973       case(ReducedFaceElements):       case(ReducedFaceElements):
974  cout << "A\n";          if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
         if (target.getFunctionSpace().getTypeCode()==FaceElements) {  
 cout << "B\n";  
975             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
 cout << "C\n";  
         throw FinleyAdapterException("A");  
976          } else {          } else {
977             throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");             throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
978         }         }
 cout << "A\n";  
979         break;         break;
980       case(Points):       case(Points):
981          if (target.getFunctionSpace().getTypeCode()==Points) {          if (target.getFunctionSpace().getTypeCode()==Points) {
# Line 1552  SystemMatrixAdapter MeshAdapter::newSyst Line 1598  SystemMatrixAdapter MeshAdapter::newSyst
1598      Paso_SystemMatrixPattern_free(fsystemMatrixPattern);      Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1599      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1600  }  }
1601    // creates a TransportProblemAdapter
1602    TransportProblemAdapter MeshAdapter::newTransportProblem(
1603                          const double theta,
1604                          const double dt_max,
1605                          const int blocksize,
1606                          const escript::FunctionSpace& functionspace,
1607                          const int type) const
1608    {
1609        int reduceOrder=0;
1610        // is the domain right?
1611        const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(functionspace.getDomain());
1612        if (domain!=*this)
1613              throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1614        // is the function space type right
1615        if (functionspace.getTypeCode()==DegreesOfFreedom) {
1616            reduceOrder=0;
1617        } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1618            reduceOrder=1;
1619        } else {
1620            throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1621        }
1622        // generate matrix:
1623        
1624        Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1625        checkFinleyError();
1626        Paso_FCTransportProblem* transportProblem;
1627        transportProblem=Paso_FCTransportProblem_alloc(theta,dt_max,fsystemMatrixPattern,blocksize);
1628        checkPasoError();
1629        Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1630        return TransportProblemAdapter(transportProblem,theta,dt_max,blocksize,functionspace);
1631    }
1632    
1633  //  //
1634  // vtkObject MeshAdapter::createVtkObject() const  // vtkObject MeshAdapter::createVtkObject() const

Legend:
Removed from v.1361  
changed lines
  Added in v.1375

  ViewVC Help
Powered by ViewVC 1.1.26