diff --git a/src/duct/duct.cpp b/src/duct/duct.cpp index 8c08eb9..ee4884b 100644 --- a/src/duct/duct.cpp +++ b/src/duct/duct.cpp @@ -231,7 +231,6 @@ vd Duct::initialSolution() const { us Duct::getNEqs() const { TRACE(15,"Duct::getNEqs()"); - us Ns = sys.Ns(); // The number of equations per gridpoint. We have: continuity, // momentum, energy, and state @@ -240,15 +239,15 @@ us Duct::getNEqs() const { // When we have to solve a solid heat balance number_eqs+= (_ductpb.stempmodel() == pb::HeatBalance ? 1 : 0); - us neqs = Ns*number_eqs*(ngp()-1); + us neqs = number_eqs*(ngp()-1); // For the last gridpoint, we also have an equation of state - neqs += Ns; + neqs += 1; // We also have an extra equation for isentropic. For the energy // transport equation, this would result in a boundary condition if(_ductpb.htmodel() == pb::Isentropic) { - neqs += Ns; + neqs += 1; } VARTRACE(15,neqs); @@ -256,14 +255,13 @@ us Duct::getNEqs() const { } us Duct::getNDofs() const { TRACE(15,"Duct::getNDofs()"); - us Ns = sys.Ns(); - + // rho,u,T,p us nvars_per_gp = 4; // Ts maybe nvars_per_gp += (_ductpb.stempmodel() == pb::HeatBalance ? 1 : 0); - return Ns*nvars_per_gp*ngp(); + return nvars_per_gp*ngp(); } d Duct::getMass() const { diff --git a/src/duct/duct.h b/src/duct/duct.h index fa0ef43..82b7ddf 100644 --- a/src/duct/duct.h +++ b/src/duct/duct.h @@ -47,7 +47,7 @@ public: * * @return The index of the degree of freedom */ - Dof getDof(int varnr,int gp) const; + PosId getDof(int varnr,int gp) const; /// Postprocessing functions vd rhox(int t) const { return getvarx(constants::rho,t);} diff --git a/src/ductbc/adiabaticwall.cpp b/src/ductbc/adiabaticwall.cpp index 10c5b88..b49d2e1 100644 --- a/src/ductbc/adiabaticwall.cpp +++ b/src/ductbc/adiabaticwall.cpp @@ -43,12 +43,12 @@ void AdiabaticWall::residualJac(arma::subview_col&& residual ) const { TRACE(15,"AdiabaticWall::residual()"); - const pb::Duct& dpb = getDuct(sys).getDuctPb(); + const Duct& duct = getDuct(); + const pb::Duct& dpb = duct.getDuctPb(); us Ns = sys.Ns(); - const Duct& duct = getDuct(sys); if(_side == pb::left) { - residual.subvec(0,Ns-1) = duct.ut(sys,0); + residual.subvec(0,Ns-1) = duct.ut(0); if(dpb.htmodel() != pb::Isentropic) { // TODO: Put the boundary condition of zero heat flux here // residual.subvec(Ns,2*Ns-1) = @@ -56,7 +56,7 @@ void AdiabaticWall::residualJac(arma::subview_col&& residual } } else { - residual.subvec(0,Ns-1) = duct.ut(sys,-1); + residual.subvec(0,Ns-1) = duct.ut(-1); if(dpb.htmodel() != pb::Isentropic) { // TODO: Put the boundary condition of zero heat flux here // residual.subvec(Ns,2*Ns-1) = duct.Tt(sys,-1) - _T; @@ -72,12 +72,12 @@ us AdiabaticWall::getNEqs() const { // dT/dx = 0 --> if htmodel is not Isentropic // dTs/dx = 0 => if stempmodel is not Prescribed - const pb::Duct& dpb = getDuct(sys).getDuctPb(); + const pb::Duct& dpb = getDuct().getDuctPb(); bool has_solideq = dpb.stempmodel() != pb::Prescribed; - us neqs = sys.Ns()*(has_solideq ? 2: 1); - if(dpb.htmodel() != pb::Isentropic) neqs+= sys.Ns(); + us neqs = (has_solideq ? 2: 1); + if(dpb.htmodel() != pb::Isentropic) neqs+= 1; VARTRACE(15,neqs); return neqs; diff --git a/src/ductbc/pressurebc.cpp b/src/ductbc/pressurebc.cpp index 0ea7ab8..25db0e6 100644 --- a/src/ductbc/pressurebc.cpp +++ b/src/ductbc/pressurebc.cpp @@ -110,8 +110,8 @@ us PressureBc::getNEqs() const { bool has_solideq = dpb.stempmodel() != pb::Prescribed; - us neqs = sys.Ns()*(has_solideq ? 2: 1); - if(dpb.htmodel() != pb::Isentropic) neqs+= sys.Ns(); + us neqs = (has_solideq ? 2: 1); + if(dpb.htmodel() != pb::Isentropic) neqs+= 1; VARTRACE(15,neqs); return neqs; diff --git a/src/sys/jaccol.cpp b/src/sys/jaccol.cpp index 6202b88..e7f0af3 100644 --- a/src/sys/jaccol.cpp +++ b/src/sys/jaccol.cpp @@ -7,63 +7,6 @@ ////////////////////////////////////////////////////////////////////// #include "jaccol.h" -#include "tasmet_variable.h" - - -JacCol::JacCol(const Variable& thevar): - coldof_(thevar.getDofNr()), - data_(thevar.getGc()->Ns(),thevar.getGc()->Ns(),fillwith::zeros) -{ } -JacCol::JacCol(us coldof,const GlobalConf& gc): - coldof_(coldof), - data_(gc.Ns(),gc.Ns(),fillwith::zeros) -{ } -JacCol::JacCol(us coldof,const dmat& data): - coldof_(coldof), - data_(data) -{ } - -JacCol::JacCol(const Variable& thevar,const dmat& data): - coldof_(thevar.getDofNr()), - data_(data) -{ } -// JacCol::JacCol(JacCol&& j): -// tobeadded(std::move(j.tobeadded)), -// coldof_(j.coldof_), -// data_(std::move(j.data_)) -// { -// TRACE(45,"JacCol::JacCol(JacCol&&)"); -// } -JacCol::JacCol(const JacCol& j): - tobeadded(j.tobeadded), - coldof_(j.coldof_), - data_(j.data_) -{ - TRACE(5,"JacCol::JacCol(const JacCol& j)"); -} -JacCol& JacCol::operator=(const JacCol& j){ - TRACE(10,"JacCol::operator=(const JacCol& j)"); - data_=j.data_; - coldof_=j.coldof_; - tobeadded=j.tobeadded; - return *this; -} -// JacCol& JacCol::operator=(JacCol&& j){ -// TRACE(45,"JacCol::operator=(const JacCol& j)"); -// data_=std::move(j.data_); -// coldof_=j.coldof_; -// tobeadded=j.tobeadded; -// return *this; -// } - -JacCol& JacCol::operator+=(const dmat& data){ - data_+=data; - return *this; -} -JacCol& JacCol::operator*=(const d& val){ - data_*=val; - return *this; -} diff --git a/src/sys/jaccol.h b/src/sys/jaccol.h index a2a287f..a963a2e 100644 --- a/src/sys/jaccol.h +++ b/src/sys/jaccol.h @@ -10,36 +10,39 @@ #define JACCOL_H #include "tasmet_types.h" -class Variable; -class GlobalConf; - +/** + * A PosId is a specifier for the position in the global Jacobian + * matrix. Both considering the column as well as the row + * + */ +struct PosId { + int segno; + int segoffset; +}; class JacCol{ // Column block of Jacobian - bool tobeadded=true; - us coldof_; // First dof of column - dmat data_; // data -public: - JacCol(const Variable&); - JacCol(us coldof,const GlobalConf&); - JacCol(us coldof,const dmat&); - JacCol(const Variable&,const dmat&); - // JacCol(JacCol&&); // Move data - JacCol(const JacCol&); - void prePostMultiply(const dmat& pre,const dmat& post){ data_=pre*data_*post;} - JacCol& operator=(const JacCol&); - // JacCol& operator=(JacCol&&); - // Negate all terms - JacCol operator-(); - bool isToAdd() const {return tobeadded;} - void setToAdd(bool set){tobeadded=set;} // If this is set to - // false, this column will not be added to the row. - JacCol& operator+=(const dmat& add); - JacCol& operator*=(const double&); - dmat& data(){return data_;} - const dmat& const_data() const {return data_;} - const us& getColDof() const{return coldof_;} - void show() const; +public: + + // For more ease of use, we make the members public. + + PosId id; /**< The Id of this DOF */ + dmat data; /**< The matrix block */ + + JacCol(PosId id): id(id){} + JacCol(PosId id,const dmat& data):id(id),data(data){} + + JacCol& operator=(const JacCol& o) { id = o.id; data=o.data; return *this; } + JacCol(const JacCol& o): id(o.id),data(o.data){} + + void prePostMultiply(const dmat& pre,const dmat& post){ data=pre*data*post;} + + // Some operations on the data + JacCol operator-() { return JacCol(id,-data);} + + JacCol& operator+=(const dmat& add) {data+=add; return *this;} + JacCol& operator*=(const double d) {data*=d; return *this;} + }; diff --git a/src/sys/jacobian.cpp b/src/sys/jacobian.cpp index 0224a54..45c3511 100644 --- a/src/sys/jacobian.cpp +++ b/src/sys/jacobian.cpp @@ -10,32 +10,34 @@ void Jacobian::operator+=(const JacRow& other){ TRACE(2,"Jacobian::append()"); jacrows.push_back(other); } -Jacobian::operator TripletList() const{ - TRACE(18,"Jacobian::operator Tripletlist()"); - int insertrow,insertcol; - TripletList res(ndofs_); - const dmat& typicaldatacel=jacrows.at(0).jaccols.at(0).const_data(); - us size=typicaldatacel.n_rows; +// Jacobian::operator TripletList() const{ - us i,j; +// TRACE(18,"Jacobian::operator Tripletlist()"); - for(const JacRow& row: jacrows) { - insertrow=row.getRowDof(); - for(const JacCol& col: row.jaccols){ - insertcol=col.getColDof(); - if(insertcol>=0){ - const dmat& data=col.const_data(); - for(i=0;i0 - } // for loop over cols - } // for loop over rows - // TRACE(10,"SFSG"); - return res; -} +// int insertrow,insertcol; +// TripletList res(ndofs_); +// const dmat& typicaldatacel=jacrows.at(0).jaccols.at(0).data; +// us size=typicaldatacel.n_rows; + +// us i,j; + +// for(const JacRow& row: jacrows) { +// insertrow=row.getRowDof(); +// for(const JacCol& col: row.jaccols){ +// insertcol=col.getColDof(); +// if(insertcol>=0){ +// const dmat& data=col.data; +// for(i=0;i0 +// } // for loop over cols +// } // for loop over rows +// // TRACE(10,"SFSG"); +// return res; +// } diff --git a/src/sys/jacobian.h b/src/sys/jacobian.h index 0c979f2..6270e84 100644 --- a/src/sys/jacobian.h +++ b/src/sys/jacobian.h @@ -14,12 +14,14 @@ class TripletList; class Jacobian{ - us ndofs_; + us _ndofs; public: - Jacobian(us ndofs):ndofs_(ndofs){} + Jacobian(us ndofs):_ndofs(ndofs){} vector jacrows; + void operator+=(const Jacobian&); void operator+=(const JacRow&); + operator TripletList() const; }; diff --git a/src/sys/jacrow.cpp b/src/sys/jacrow.cpp index a5bcadd..3b89247 100644 --- a/src/sys/jacrow.cpp +++ b/src/sys/jacrow.cpp @@ -3,48 +3,30 @@ #include "globalconf.h" - -JacRow::JacRow(const JacCol& j): - JacRow(-1,1) -{ - (*this)+=j; -} -JacRow::JacRow(us rowdof,const JacCol& j): - JacRow(rowdof,1) -{ - (*this)+=j; -} -// JacRow& JacRow::operator+=(JacCol&& j){ -// TRACE(45,"JacRow::operator+=(JacCol&& j)"); -// jaccols.emplace_back(std::move(j)); -// } JacRow& JacRow::operator+=(const JacCol& j){ TRACE(10,"JacRow::operator+=(const JacCol& j)"); - if(j.isToAdd()) - jaccols.emplace_back(j); + jaccols.push_back(j); return *this; } - -// JacRow& JacRow::addCol(const JacCol& jaccol){ -// if(jaccol.isToAdd()) -// jaccols.push_back(jaccol); -// return *this; -// } -JacRow& JacRow::operator*=(const d& val){ + +JacRow& JacRow::operator*=(const d val){ TRACE(2,"Jacobian::operator*=()"); for(JacCol& col: jaccols) - col.data()*=val; + col*=val; return *this; } JacRow& JacRow::operator+=(const JacRow& other){ - TRACE(2,"Jacobian::operator*=()"); + // Reserve some space this->jaccols.reserve(this->jaccols.capacity()+other.jaccols.size()-this->jaccols.size()); + for(const JacCol& col :other.jaccols) (*this)+=col; + return *this; } JacRow JacRow::operator-() const{ TRACE(15,"JacRow::operator-()"); + JacRow result(*this); for (JacCol& jaccol : result.jaccols) jaccol*=-1; diff --git a/src/sys/jacrow.h b/src/sys/jacrow.h index bad5387..2716956 100644 --- a/src/sys/jacrow.h +++ b/src/sys/jacrow.h @@ -8,32 +8,43 @@ #pragma once #ifndef JACROW_H #define JACROW_H - +#include #include "jaccol.h" -class Variable; - - class JacRow{ // Row in Jacobian matrix - int rowdof_=-1; // Number of first row, default is - // invalid state + public: - // Negate all terms + + PosId id; /**< Specifies the row for this + JacRow */ + std::vector jaccols; /**< The columns in this row */ + + /** + * Initialize a Jacobian row with a single column block + * + * @param id The equation ID + * @param col The added column + */ + JacRow(PosId id,const JacCol& col):JacRow(id){jaccols.push_back(col);} + + JacRow(PosId id): id(id){} + JacRow operator-() const; - vector jaccols; // Column blocks - JacRow(const JacCol&); - JacRow(us rowdof,const JacCol&); // Initialization with only one column - JacRow(int rowdofnr,us cols=2): rowdof_(rowdofnr){ jaccols.reserve(cols);} - // void addCol(const JacCol& jaccol); - // JacRow& operator+=(JacCol&&); JacRow& operator+=(const JacCol&); JacRow& operator+=(const JacRow& jacrow); - JacRow& operator*=(const d& val); // Multiply all terms with constant value + JacRow& operator*=(const d val); // Multiply all terms with constant value + + /** + * Pre- and postmultiply all Column blocks with the given + * matrices. Hence the new columns will be col_new = pre*col_old*post + * + * @param pre : The pre-multiply matrix + * @param post: The post-multiply matrix + */ void prePostMultiply(const dmat& pre,const dmat& post); - const int& getRowDof() const {return rowdof_;} - void setRowDof(us dofnr){rowdof_=dofnr;} - void show() const; + + };