Select Git revision
gsl_newt.cpp
Forked from
gmsh / gmsh
Source project has a limited visibility.
-
Jean-François Remacle authoredJean-François Remacle authored
solverAlgorithms.h 9.50 KiB
// Gmsh - Copyright (C) 1997-2013 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
//
// Contributor(s):
// Eric Bechet
//
#ifndef _SOLVERALGORITHMS_H_
#define _SOLVERALGORITHMS_H_
#include "dofManager.h"
#include "terms.h"
#include "quadratureRules.h"
#include "MVertex.h"
template<class Iterator, class Assembler> void Assemble(BilinearTermBase &term, FunctionSpaceBase &space,
Iterator itbegin, Iterator itend,
QuadratureBase &integrator, Assembler &assembler)
// symmetric
{
fullMatrix<typename Assembler::dataMat> localMatrix;
std::vector<Dof> R;
for (Iterator it = itbegin; it != itend; ++it){
MElement *e = *it;
R.clear();
IntPt *GP;
int npts = integrator.getIntPoints(e, &GP);
term.get(e, npts, GP, localMatrix); //localMatrix.print();
space.getKeys(e, R);
assembler.assemble(R, localMatrix);
}
}
template<class Assembler> void Assemble(BilinearTermBase &term, FunctionSpaceBase &space, MElement *e,
QuadratureBase &integrator, Assembler &assembler) // symmetric
{
fullMatrix<typename Assembler::dataMat> localMatrix;
std::vector<Dof> R;
IntPt *GP;
int npts = integrator.getIntPoints(e, &GP);
term.get(e, npts, GP, localMatrix);
space.getKeys(e, R);
assembler.assemble(R, localMatrix);
}
template<class Iterator, class Assembler> void Assemble(BilinearTermBase &term,
FunctionSpaceBase &shapeFcts,
FunctionSpaceBase &testFcts,
Iterator itbegin, Iterator itend,
QuadratureBase &integrator,
Assembler &assembler) // non symmetric
{
fullMatrix<typename Assembler::dataMat> localMatrix;
std::vector<Dof> R;
std::vector<Dof> C;
for (Iterator it = itbegin; it != itend; ++it){
MElement *e = *it;
R.clear();
C.clear();
IntPt *GP;
int npts = integrator.getIntPoints(e, &GP);
term.get(e, npts, GP, localMatrix); //localMatrix.print();
printf("local matrix size = %d %d\n", localMatrix.size1(), localMatrix.size2());
shapeFcts.getKeys(e, R);
testFcts.getKeys(e, C);// std::cout << "assembling normal test function ; lagrange trial function : " << std::endl;
// for (int i = 0 ; i < R.size() ; i++)
// {
// std::cout << "tests : " << R[i].getEntity() << ":" << R[i].getType() << std::endl ;
// }
// for (int i = 0 ; i < R.size() ; i++)
// {
// std::cout << "trial : " << C[i].getEntity() << ":" << C[i].getType() << std::endl ;
// }
assembler.assemble(R, C, localMatrix);
// std::cout << "assembling lagrange test function ; normal trial function : " << std::endl;
// for (int i = 0 ; i < R.size() ; i++)
// {
// std::cout << "tests : " << C[i].getEntity() << ":" << C[i].getType() << std::endl ;
// }
// for (int i = 0 ; i < R.size() ; i++)
// {
// std::cout << "trial : " << R[i].getEntity() << ":" << R[i].getType() << std::endl ;
// }
assembler.assemble(C, R, localMatrix.transpose());
}
}
template<class Iterator, class Assembler> void Assemble(LinearTermBase<double> &term, FunctionSpaceBase &space,
Iterator itbegin, Iterator itend,
QuadratureBase &integrator, Assembler &assembler)
{
fullVector<typename Assembler::dataMat> localVector;
std::vector<Dof> R;
for (Iterator it = itbegin; it != itend; ++it){
MElement *e = *it;
R.clear();
IntPt *GP;
int npts = integrator.getIntPoints(e, &GP);
term.get(e, npts, GP, localVector); //localVector.print();
space.getKeys(e, R);
assembler.assemble(R, localVector);
}
}
template<class Assembler> void Assemble(LinearTermBase<double> &term, FunctionSpaceBase &space, MElement *e,
QuadratureBase &integrator, Assembler &assembler)
{
fullVector<typename Assembler::dataMat> localVector;
std::vector<Dof> R;
IntPt *GP;
int npts = integrator.getIntPoints(e, &GP);
term.get(e, npts, GP, localVector);
space.getKeys(e, R);
assembler.assemble(R, localVector);
}
template<class Iterator, class dataMat> void Assemble(ScalarTermBase<double> &term,
Iterator itbegin, Iterator itend,
QuadratureBase &integrator, dataMat & val)
{
dataMat localval;
for (Iterator it = itbegin; it != itend; ++it){
MElement *e = *it;
IntPt *GP;
int npts = integrator.getIntPoints(e, &GP);
term.get(e, npts, GP, localval);
val += localval;
}
}
template<class Iterator, class dataMat> void Assemble(ScalarTermBase<double> &term, MElement *e,
QuadratureBase &integrator, dataMat & val){
dataMat localval;
IntPt *GP;
int npts = integrator.getIntPoints(e, &GP);
term.get(e, npts, GP, localval);
val += localval;
}
template<class Assembler> void FixDofs(Assembler &assembler, std::vector<Dof> &dofs,
std::vector<typename Assembler::dataVec> &vals)
{
int nbff = dofs.size();
for (int i = 0; i < nbff; ++i){
assembler.fixDof(dofs[i], vals[i]);
}
}
class FilterDof
{
public:
virtual ~FilterDof(){}
virtual bool operator()(Dof key) = 0;
};
class FilterDofTrivial : public FilterDof
{
public:
virtual bool operator()(Dof key) {return true;}
};
class FilterDofComponent : public FilterDof
{
int comp;
public :
FilterDofComponent(int comp_) : comp(comp_) {}
virtual bool operator()(Dof key)
{
int type = key.getType();
int icomp, iphys;
Dof::getTwoIntsFromType(type, icomp, iphys);
if (icomp == comp) return true;
return false;
}
};
// return true if the Dof is in a filled set
class FilterDofSet : public FilterDof
{
protected:
std::set<Dof> _dofset;
public:
FilterDofSet() : FilterDof(){}
virtual ~FilterDofSet(){}
virtual bool operator()(Dof key)
{
std::set<Dof>::iterator itR = _dofset.find(key);
if(itR == _dofset.end()){
return false;
}
return true;
}
virtual void addDof(Dof key)
{
_dofset.insert(key);
}
virtual void addDof(std::vector<Dof> &R)
{
for(unsigned int i=0;i<R.size();i++)
this->addDof(R[i]); }
};
template<class Assembler> void FixNodalDofs(FunctionSpaceBase &space, MElement *e, Assembler &assembler,
simpleFunction<typename Assembler::dataVec> &fct,
FilterDof &filter)
{
std::vector<MVertex*> tabV;
int nv = e->getNumVertices();
std::vector<Dof> R;
space.getKeys(e, R);
tabV.reserve(nv);
for (int i = 0; i < nv; ++i)
tabV.push_back(e->getVertex(i));
for (std::vector<Dof>::iterator itd = R.begin(); itd != R.end(); ++itd){
Dof key = *itd;
if (filter(key)){
for (int i = 0;i < nv; ++i){
if (tabV[i]->getNum() == key.getEntity()){
assembler.fixDof(key, fct(tabV[i]->x(), tabV[i]->y(), tabV[i]->z()));
break;
}
}
}
}
}
template<class Iterator, class Assembler> void FixNodalDofs(FunctionSpaceBase &space, Iterator itbegin,
Iterator itend, Assembler &assembler,
simpleFunction<typename Assembler::dataVec> &fct,
FilterDof &filter)
{
for (Iterator it = itbegin; it != itend; ++it)
FixNodalDofs(space, *it, assembler, fct, filter);
}
template<class Iterator, class Assembler> void FixVoidNodalDofs(FunctionSpaceBase &space, Iterator itbegin,
Iterator itend, Assembler &assembler)
{
FilterDofTrivial filter;
simpleFunction<typename Assembler::dataVec> fct(0.);
FixNodalDofs(space, itbegin, itend, assembler, fct, filter);
}
template<class Iterator, class Assembler> void NumberDofs(FunctionSpaceBase &space, Iterator itbegin,
Iterator itend, Assembler &assembler)
{
for (Iterator it = itbegin; it != itend; ++it){
MElement *e = *it;
std::vector<Dof> R;
space.getKeys(e, R);
int nbdofs = R.size();
for (int i = 0; i < nbdofs; ++i) assembler.numberDof(R[i]);
}
}
//// Mean HangingNodes
//template <class Assembler> void FillHangingNodes(FunctionSpaceBase &space, std::map<int,std::vector <int> > &HangingNodes, Assembler &assembler, int &field, int &dim)
//{
// std::map<int, std::vector <int> >::iterator ith;
// ith = HangingNodes.begin();
// int compt = 1;
// while (ith != HangingNodes.end()){
// float fac;
// fac = 1.0 / (ith->second).size();
// for (int j = 0; j < dim; j++){
// DofAffineConstraint<double> constraint;
// int type = Dof::createTypeWithTwoInts(j, field);
// Dof hgnd(ith->first, type);// for (int i = 0; i < (ith->second).size(); i++){
// Dof key((ith->second)[i], type);
// std::pair<Dof, double > linDof(key, fac);
// constraint.linear.push_back(linDof);
// }
// constraint.shift = 0;
// assembler.setLinearConstraint (hgnd, constraint);
// }
// ith++;
// compt++;
// }
//}
#endif// _SOLVERALGORITHMS_H_