библиотека nla3d с улучшенным выводом VTK-файлов (3 этап)

README.md

@@ -6,4 +6,8 @@
 2. С помощью проекта CodeBlocks (файл math.cbp) или make-файла откомпилировать библиотеку math, находящуюся в папке math.
 По умолчанию включена поддержка как решателей paralution (-DNLA3D_USE_PARALUTION), так и решателей mkl (-DNLA3D_USE_MKL). При необходимости, ненужное отключить в соответствующих настройках проекта CodeBlocks или в make-файле.
 Особое внимание уделить файлу EquationSolver.h и значениям абсолютной невязки "double abs_tol = 1.0E-06" и относительной невязки "double rel_tol = 1.0E-06" для решателей paralution. Это критерии выхода для итерационного решателя. В зависимости от требуемой точности решения их нужно изменять.
+3. С помощью проекта CodeBlocks (файл nla3d.cbp) или make-файла откомпилировать библиотеку nla3d, находящуюся в папке nla3d.
+В данной версии nla3d улучшен алгоритм вывода данных для VTK-файлов. Скорость вывода увеличена в разы, особенно хорошо это заметно на больших сетках.
+Для успешной компиляции необходимо наличие исходников библиотеки линейной алгебры eigen.
+В ходе компиляции будет большое количество warning'ов, это нормально и связано с eigen.
 

nla3d/Dof.cpp

+#include "Dof.h"
+
+namespace nla3d {
+
+const uint16 Dof::numberOfDofTypes = Dof::UNDEFINED;
+const char* const Dof::dofTypeLabels[] = {"UX", "UY", "UZ", "ROTX", "ROTY",
+  "ROTZ", "HYDRO_PRESSURE", "TEMP", "UNDEFINED"};
+
+
+// Check that TypeLabels has the same length as numberOfDofTypes
+static_assert(Dof::numberOfDofTypes == sizeof(Dof::dofTypeLabels)/sizeof(Dof::dofTypeLabels[0]) - 1,
+                "dofTypeLabels and dofType must have the same number of elements");
+
+Dof::dofType Dof::label2dofType (const std::string& label) {
+  for (uint16 i = 0; i < numberOfDofTypes; i++) {
+    if (label.compare(dofTypeLabels[i]) == 0) {
+      return (dofType) i;
+    }
+  }
+  LOG(ERROR) << "Unknown dof key = " << label;
+  return Dof::UNDEFINED;
+}
+
+
+void DofCollection::initDofTable(uint32 _numberOfEntities) {
+  clearDofTable();
+  numberOfEntities = _numberOfEntities;
+  dofPos.assign(numberOfEntities + 1, 0);
+}
+
+
+// n index starts from 1
+void DofCollection::addDof(uint32 n, std::initializer_list<Dof::dofType> __dofs) {
+  assert(n <= numberOfEntities);
+  assert(dofPos.size() > 0);
+  std::set<Dof::dofType> _dofs(__dofs);
+  std::vector<Dof> newDofs;
+  // if already registered - just return
+  if (dofPos[n-1] == dofPos[n]) {
+    for (auto v : _dofs)
+      newDofs.push_back(Dof(v));
+  } else {
+    for (auto v : _dofs) {
+      bool isFound = false;
+      for (auto it = dofs.begin() + dofPos[n-1]; it < dofs.begin() + dofPos[n]; it++) {
+        if (it->type == v) {
+          isFound = true;
+          break;
+        }
+      }
+      if (!isFound) {
+        newDofs.push_back(Dof(v));
+      }
+    }
+  }
+
+  if (newDofs.size() == 0) return;
+
+  dofs.insert(dofs.begin() + dofPos[n], newDofs.begin(), newDofs.end());
+  uniqueDofTypes.insert(std::begin(_dofs), std::end(_dofs));
+
+  // incement dofPos with newDofs.size() for all above
+  uint32 i = n;
+  while (i <= numberOfEntities) dofPos[i++] += newDofs.size();
+
+  numberOfUsedDofs += newDofs.size();
+}
+
+
+void DofCollection::clearDofTable() {
+  dofs.clear();
+  dofPos.clear();
+  uniqueDofTypes.clear();
+  numberOfUsedDofs = 0;
+  numberOfEntities = 0;
+}
+
+
+// n index starts from 1
+bool DofCollection::isDofUsed(uint32 n, Dof::dofType dof) {
+  assert(n <= numberOfEntities);
+  assert(dofPos.size() > 0);
+  if (dofPos[n-1] == dofPos[n]) return false;
+  for (auto it = dofs.begin() + dofPos[n-1]; it < dofs.begin() + dofPos[n]; it++) {
+    if (it->type == dof) return true;
+  }
+  return false;
+}
+
+
+// return nullptr if Dof was not found
+Dof* DofCollection::getDof(uint32 n, Dof::dofType dof) {
+  assert(n <= numberOfEntities);
+  assert(dofPos.size() > 0);
+
+  Dof* pdof = nullptr;
+
+  for (auto it = dofs.begin() + dofPos[n-1]; it < dofs.begin() + dofPos[n]; it++) {
+    if (it->type == dof) {
+      pdof = &(*it);
+      break;
+    }
+  }
+  return pdof;
+}
+
+
+std::set<Dof::dofType> DofCollection::getUniqueDofTypes() {
+  return uniqueDofTypes;
+}
+
+
+} // namespace nla3d

nla3d/Dof.h

+// This file is a part of nla3d project. For information about authors and
+// licensing go to project's repository on github:
+// https://github.com/dmitryikh/nla3d 
+
+#pragma once
+#include "sys.h"
+
+namespace nla3d {
+
+class DofCollection;
+// class for Degree of Freedom informations 
+// (is it fixed boundary condition, what its number in equation system, ..)
+class Dof {
+  public:
+    enum dofType {
+      UX = 0,
+      UY,
+      UZ,
+      ROTX,
+      ROTY,
+      ROTZ,
+      HYDRO_PRESSURE,
+      TEMP,
+      UNDEFINED // should be last
+    };
+
+    Dof(dofType t) : eqNumber(0), isConstrained(false), type(t) { }
+
+    uint32 eqNumber = 0; // 0 - not use
+    bool isConstrained = false; // is this DOF defined by B.C.
+    dofType type;
+
+    static const char* const dofTypeLabels[];
+    static const uint16 numberOfDofTypes;
+    static dofType label2dofType (const std::string& label);
+    static const char* dofType2label (const dofType t);
+};
+
+inline const char* Dof::dofType2label(const Dof::dofType t) {
+  assert(t < Dof::UNDEFINED);
+  return dofTypeLabels[t];
+}
+
+
+class DofCollection {
+  public:
+    uint32 getNumberOfUsedDofs();
+    uint32 getNumberOfEntities();
+
+    // return a poiner to Dof class for Entity n and for type dof
+    Dof* getDof(uint32 n, Dof::dofType dof);
+
+    // return begin and end iterator of Dofs for Entity n
+    std::pair<std::vector<Dof>::iterator,
+              std::vector<Dof>::iterator> getEntityDofs(uint32 n);
+
+    void initDofTable(uint32 _numberOfEntities);
+    void addDof(uint32 n, std::initializer_list<Dof::dofType> _dofs);
+    void clearDofTable();
+    bool isDofUsed(uint32 n, Dof::dofType dof);
+
+    std::set<Dof::dofType> getUniqueDofTypes();
+
+  private:
+    uint32 numberOfUsedDofs = 0;
+    uint32 numberOfEntities = 0;
+
+    // vector of Dof objects 
+    std::vector<Dof> dofs;
+    // array of indexes to find where dofs for particular entity is located in dofs
+    // Dof for entity n will be located from dofPos[n-1] included to dofPos[n] excluded 
+    std::vector<uint32> dofPos;
+    // set of unique dofs used in collection
+    std::set<Dof::dofType> uniqueDofTypes;
+};
+
+
+inline uint32 DofCollection::getNumberOfUsedDofs() {
+  return numberOfUsedDofs;
+}
+
+
+inline uint32 DofCollection::getNumberOfEntities() {
+  return numberOfEntities;
+}
+
+
+inline std::pair<std::vector<Dof>::iterator,
+                 std::vector<Dof>::iterator> DofCollection::getEntityDofs(uint32 n) {
+  assert(n > 0);
+  assert(n <= numberOfEntities);
+  assert(dofPos.size() > 0);
+
+  return std::make_pair(dofs.begin() + dofPos[n-1], dofs.begin() + dofPos[n]);
+}
+
+
+} // namespace nla3d 

nla3d/FEComponent.cpp

+// This file is a part of nla3d project. For information about authors and
+// licensing go to project's repository on github:
+// https://github.com/dmitryikh/nla3d 
+
+#include "FEComponent.h"
+
+namespace nla3d {
+
+const char* const FEComponent::labelsOfComponent[]={"NOT DEFINED","ELEM", "NODE", "LAST"};
+
+FEComponent::FEComponent () {
+  type = NOT_DEFINED;
+  name = "";
+}
+FEComponent::~FEComponent () {
+  list.clear();
+} 
+
+FEComponent::typeOfComponent FEComponent::typeFromString(const std::string& typeName) {
+  for (size_t i = 1; i < LAST; i++) {
+    if (typeName.compare(labelsOfComponent[i]) == 0) {
+      return static_cast<typeOfComponent> (i);
+    }
+  }
+  LOG(ERROR) << "Can't find FEComponent::type by name " << typeName;
+  return NOT_DEFINED;
+}
+
+} // namespace nla3d

nla3d/FEComponent.h

+// This file is a part of nla3d project. For information about authors and
+// licensing go to project's repository on github:
+// https://github.com/dmitryikh/nla3d 
+
+#pragma once
+#include "sys.h"
+
+namespace nla3d {
+
+class FEComponent {
+  public:
+    enum typeOfComponent {
+      NOT_DEFINED,
+      ELEMENTS,
+      NODES,
+      LAST
+    };
+
+    static const char* const labelsOfComponent[];
+
+    FEComponent ();
+    ~FEComponent ();
+
+    std::vector<uint32> list;
+    typeOfComponent type;
+    std::string name;
+
+    static typeOfComponent typeFromString(const std::string& typeName);
+    void print ();
+};
+
+inline MAKE_LOGGABLE(FEComponent, obj, os) {
+  os << obj.name << ": " << obj.list.size() << " " << obj.labelsOfComponent[obj.type];
+  return os;
+} 
+
+} // namespace nla3d 

nla3d/FESolver.h

+// This file is a part of nla3d project. For information about authors and
+// licensing go to project's repository on github:
+// https://github.com/dmitryikh/nla3d 
+
+#pragma once
+#include "sys.h"
+#include "math/Vec.h"
+#include "math/EquationSolver.h"
+#include "FEStorage.h"
+#include "PostProcessor.h"
+
+#include <Eigen/Dense>
+
+namespace nla3d {
+
+using namespace math;
+
+class FEStorage;
+class PostProcessor;
+
+// TimeControl class helps to control time stepping and equilibrium iterations.
+// Before solution one can setup setStartTime() and setEndTime() and then perform time stepping with
+// nextStep(delta) where delta is solver specific time step. For every time step many equilibrium
+// iterations could be performed by nextEquilibriumStep(). All intermediate time steps are stored in
+// convergedTimeInstances, for every time steps number of equilibrium steps are stored in
+// equilibriumSteps.
+class TimeControl {
+public:
+  uint16 getCurrentStep ();
+  uint16 getNumberOfConvergedSteps ();
+
+  uint16 getCurrentEquilibriumStep ();
+  uint16 getTotalNumberOfEquilibriumSteps ();
+
+  bool nextStep (double delta);
+  void nextEquilibriumStep ();
+
+  double getCurrentTime ();
+  double getCurrentNormalizedTime ();
+  double getEndTime ();
+  double getStartTime ();
+  double getCurrentTimeDelta ();
+  double getCurrentNormalizedTimeDelta ();
+
+  void setEndTime (double _endTime);
+  void setStartTime (double _startTime);
+protected:
+  std::list<double> convergedTimeInstances;
+  std::list<uint16> equilibriumSteps;
+  uint16 currentEquilibriumStep = 0;
+  uint16 totalNumberOfEquilibriumSteps = 0;
+  double currentTimeDelta = 0.0;
+  double endTime = 1.0;
+  double startTime = 0.0;
+  double currentTime = 0.0;
+};
+
+// The abstract class that represents the FE solver. This class use information and methods from
+// FEStorage in order to apply particular solution scheme. For example, here are already some of
+// them: solver for linear steady FE model with concentrated loads and fixed DoFs (LinearFESolver),
+// solver for non-linear quasi-steady FE model based on incremental approach (NonlinearFESolver),
+// solver for linear transient FE model based on Newmark time-integration scheme
+// (LinearTransientFESolver).
+//
+// FESolver keeps a pointer to FEStorage to receive assembled global equation system, and a pointer
+// to math::EquationSolver to solve linear system of equations.
+//
+// FESovler manages PostProcessors objects in order to call them on every time step to do some
+// useful routines (write res files, log reaction forces, etc.)
+//
+// FESolver by default contains arrays of loadBC and fixBC structures in order to apply constant BC
+// on DoFs (concentrated load and fixed DoF value). But particular realizations of FESolver can
+// incorporate time-depended BC along with or instead of loadBC and fixBC.
+//
+// FESolver by mean of initSolutionData() keeps pointers to FEStorage matrices and vectors. matK for
+// stiffness matrix, matC for dumping matrix, matM for inertia matrix. vecU - for DoF values vector
+// with references to particular parts of it: vecUc(part for constrained DoFs), vecUs(part for
+// unknown DoFs), vecUl(part for Lagrangian lambdas from Mpc equations), vecUsl(combination of vecUs and
+// vecUl). The same for first and second time derivatives: vecDU, vecDDU.
+// Also there are references to RHS parts: vecF - for FE internal loads. vecFl has special treatment
+// as RHS of Mpc equations. vecR - for external concentrated loads. vecRc is treated as reaction
+// forces of fixation for constrained DoFs, vecRs is treated as external concetrated loads and vecRl
+// should be zero all times. All references above are just easy way to access FEStorage internal
+// data. Actually all this vectors and matrices ale allocated inside of FEStorage. By performing
+// FEStorage::assembleGlobalEqMatrices() matK/C/M, vecF are assembled with particular numbers based
+// on FE element formulation and on Mpc equations. Other entities (vecU, vecDU, vecDDU, vecR) are
+// fully under FESovler control. FESolver is responsible on timely updating this values.
+//
+class FESolver {
+  public:
+    FESolver ();
+    virtual ~FESolver ();
+
+    void attachFEStorage(FEStorage *st);
+    void attachEquationSolver(math::EquationSolver *eq);
+
+    // main method were all solution scheme specific routines are performed
+    virtual void solve() = 0;
+
+    size_t getNumberOfPostProcessors();
+    // get an instance of PostProcessor by its number
+    // _np >= 0
+    PostProcessor& getPostProcessor(size_t _np);
+    // The function stores PostProcessor pointer in FESolver. FESovler will free the memory by itself.
+    uint16 addPostProcessor(PostProcessor *pp);
+    void deletePostProcessors();
+
+    // run FEStorage procedures for initialization solution data structures and map matK, matC, matM,
+    // vecU, vecDU, vecDDU, vecR, vecF pointer on FEStorage's allocated structures.
+    void initSolutionData();
+
+    // The function applies boundary conditions stored in `loads` and `fixs` arrays.
+    // `time` should be normalized(in range [0.0; 1.0])
+    void applyBoundaryConditions(double time);
+
+    // constrained DoFs has special treatment in FEStorage, that's why before initialization of
+    // solution data we need to tell FEStorage which DoFs is constrained. setConstrainedDofs() does
+    // this work based on `fixs` array.
+    void setConstrainedDofs();
+
+    // add DoF fixation (constraint) boundary condition 
+    void addFix(int32 n, Dof::dofType dof, const double value = 0.0);
+    // add DoF load (force) boundary condition 
+    void addLoad(int32 n, Dof::dofType dof, const double value = 0.0);
+    
+    // for debug purpose:
+    // dump matrices matK, matC, matM and vectors vecF, vecR
+    void dumpMatricesAndVectors(std::string filename);
+    // read matrices matK, matC, matM and vectors vecF, vecR from file and compare them with
+    // solution instances
+    void compareMatricesAndVectors(std::string filename, double th = 1.0e-9);
+  protected:
+    FEStorage* storage = nullptr;
+    math::EquationSolver* eqSolver = nullptr;
+
+    // Array of PostProcessors. Here is a conception of PostProcessors that can do useful work every
+    // iteration of the solution. For example here is a VtkProcessor class to write *.vtk file with
+    // model and solution data (displacements, stresses and others). An instance of PostProcessor
+    // class is created outside of FESolver. But with function addPostProcessor(..) FESolver take
+    // control on the PostProcessor instance.  The memory released in deletePostProcessors().
+    std::vector<PostProcessor*> postProcessors;
+
+    // the references on FEStorage FE data structures which is frequently used by FESolver. This
+    // references make it easy to operate with important FE data structures. 
+    BlockSparseSymMatrix<2>* matK = nullptr;
+    BlockSparseSymMatrix<2>* matC = nullptr;
+    BlockSparseSymMatrix<2>* matM = nullptr;
+
+    // vecU is a reference on a whole DoF values vector, but vecUc, vecUs, vecUl, vecUsl are
+    // references on parts of the full vector. Some times it's handy to operate with the particular
+    // part only.
+    dVec vecU;
+    dVec vecUc;
+    dVec vecUs;
+    dVec vecUl;
+    dVec vecUsl;
+
+    // for first time derivatives
+    dVec vecDU;
+    dVec vecDUc;
+    dVec vecDUs;
+    dVec vecDUl;
+    dVec vecDUsl;
+
+    // for second time derivatives
+    dVec vecDDU;
+    dVec vecDDUc;
+    dVec vecDDUs;
+    dVec vecDDUl;
+    dVec vecDDUsl;
+
+    // for internal element loads
+    dVec vecF;
+    dVec vecFc;
+    dVec vecFs;
+    dVec vecFl;
+    dVec vecFsl;
+
+    // for external loads
+    dVec vecR;
+    dVec vecRc;
+    dVec vecRs;
+    dVec vecRl;
+    dVec vecRsl;
+
+    // List of concentrated load boundary conditions (addition to RHS of global eq. system)
+    std::list<loadBC> loads;
+
+    // List of prescribed values for DoFs.
+    // NOTE: DoFs with fixed values are treated in nla3d in a special manner: this DoFs are eliminated from
+    // global solve system and then reaction loads (vecRc vector) of a such DoFs are found.
+    std::list<fixBC> fixs;
+};
+
+
+// particular realization of FESolver for linear tasks
+class LinearFESolver : public FESolver {
+  public:
+    LinearFESolver();
+    virtual void solve();
+};
+
+// Iterative solver for Full Newton-Raphson procedure
+// The convergence is controlled by mean increment of DoF values
+class NonlinearFESolver : public FESolver {
+  public:
+    NonlinearFESolver();
+
+    TimeControl timeControl;
+    uint16 numberOfIterations = 20;
+    uint16 numberOfLoadsteps = 10;
+
+    double convergenceCriteria = 1.0e-3;
+
+    virtual void solve();
+  protected:
+    double calculateCriteria(dVec& delta);
+};
+
+
+// Solver for time integration of linear systems: M * DDU + C * DU + K * U = F + R
+// use Newmark scheme (based on section 9.2.4. Bathe K.J., Finite Element Procedures, 1997)
+class LinearTransientFESolver : public FESolver {
+  public:
+    LinearTransientFESolver();
+
+    uint16 numberOfTimesteps = 100;
+
+    double alpha = 0.25; // trapezoidal rule by default
+    double delta = 0.5;
+
+    // start time
+    double time0 = 0.0;
+    // end time
+    double time1 = 1.0;
+    // initial values for vecUc
+    double initValue = 0.0;
+
+    virtual void solve ();
+
+    double a0, a1, a2, a3, a4, a5, a6, a7;
+};
+
+} // namespace nla3d

nla3d/Mpc.cpp

+// This file is a part of nla3d project. For information about authors and
+// licensing go to project's repository on github:
+// https://github.com/dmitryikh/nla3d 
+
+#include "Mpc.h"
+#include "FEStorage.h"
+#include "Node.h"
+
+namespace nla3d {
+
+using namespace math;
+
+
+void Mpc::print (std::ostream& out) {
+  list<MpcTerm>::iterator token = eq.begin();
+  bool first = true;
+  while (token != eq.end()) {
+    //Cij_add(eq_num, token->node, token->node_dof, token->coef);
+    if (!first) {
+      out << " + ";
+    } else {
+      first = false;
+    }
+    out << toStr(token->coef) << " * " << toStr(token->node) << ":" << Dof::dofTypeLabels[token->node_dof];
+    token++;
+  }
+  out << " = " << toStr(b);
+}
+
+void MpcCollection::registerMpcsInStorage () {
+  for (size_t i = 0; i < collection.size(); i++) {
+    storage->addMpc(collection[i]);
+  }
+}
+
+
+void MpcCollection::printEquations (std::ostream& out) {
+  for (size_t i = 0; i < collection.size(); i++) {
+    collection[i]->print(out);
+    out << std::endl;
+  }
+}
+
+RigidBodyMpc::RigidBodyMpc () {
+
+}
+
+RigidBodyMpc::~RigidBodyMpc () {
+  collection.clear();
+}
+
+void RigidBodyMpc::pre () {
+  // create Mpc for all nodes and all dofs..
+  collection.assign(slaveNodes.size() * dofs.size(), NULL);
+  for (size_t i = 0; i < slaveNodes.size(); i++) {
+    for (size_t j = 0; j < dofs.size(); j++) {
+      Mpc* mpc = new Mpc;
+      mpc->b = 0.0;
+      mpc->eq.push_back(MpcTerm(slaveNodes[i], dofs[j], 1.0));
+      mpc->eq.push_back(MpcTerm(masterNode, dofs[j], -1.0));
+      mpc->eq.push_back(MpcTerm(masterNode, Dof::ROTX, 0.0));
+      mpc->eq.push_back(MpcTerm(masterNode, Dof::ROTY, 0.0));
+      mpc->eq.push_back(MpcTerm(masterNode, Dof::ROTZ, 0.0));
+      collection[i*dofs.size() + j] = mpc;
+    }
+  }
+  storage->addNodeDof(masterNode, {Dof::UX, Dof::UY, Dof::UZ, Dof::ROTX, Dof::ROTY, Dof::ROTZ});
+}
+
+void RigidBodyMpc::update () {
+  // rigid body motion formula is the next:
+  // {u} = {w} +([C]-[I])*({p}-{q})
+  // where  {u} - displacement of a slave node
+  //        {w} - displacement of the master node
+  //        [C] - rotation matrix (depends on rotation vector {theta})
+  //        {p} - position of the master node
+  //        {q} - position of a slave node
+  // where [C] in componentwise form:
+  // C_{ij} = \delta_{ij} cos\theta + \frac{sin\theta}{\theta}e_{ikj}\theta_k + \left(\frac{1-cos\theta}{\theta^2}\theta_i \theta_j\right) 
+  Vec<3> theta0;
+  Vec<3> w0;
+  theta0[0] = storage->getNodeDofSolution(masterNode, Dof::ROTX);
+  theta0[1] = storage->getNodeDofSolution(masterNode, Dof::ROTY);
+  theta0[2] = storage->getNodeDofSolution(masterNode, Dof::ROTZ);
+  w0[0] = storage->getNodeDofSolution(masterNode, Dof::UX);
+  w0[1] = storage->getNodeDofSolution(masterNode, Dof::UY);
+  w0[2] = storage->getNodeDofSolution(masterNode, Dof::UZ);
+  double thetaNorm = theta0.length();
+  Vec<3> masterPos;
+  Vec<3> slavePos;
+  Mat<3,3> C;
+  double c1, c2, c3, c4;
+  if (thetaNorm < 1.0e-3) {
+    // for very small angles we need to define some undefined values
+    // TODO: for now eps = 1.0e-3 was choosed randomly, need to check it
+    c1 = 1.0;
+    c2 = 0.5;
+    c3 = -1.0/3.0;
+    c4 = -1.0/12.0;
+  } else {
+  // define constant exactly
+    c1 = sin(thetaNorm) / thetaNorm;
+    c2 = (1.0 - cos(thetaNorm)) / (thetaNorm * thetaNorm);
+    c3 = (thetaNorm * cos(thetaNorm) - sin(thetaNorm)) / (thetaNorm * thetaNorm * thetaNorm);
+    c4 = (thetaNorm * sin(thetaNorm) - 2.0 + 2.0 * cos(thetaNorm)) / (thetaNorm * thetaNorm * thetaNorm * thetaNorm);
+  }
+  for (uint16 i = 0; i < 3; i++) {
+    for (uint16 j = 0; j < 3; j++) {
+      double theta0Mat_ij = solidmech::LeviCivita[i][0][j] * theta0[0] +
+                      solidmech::LeviCivita[i][1][j] * theta0[1] +
+                      solidmech::LeviCivita[i][2][j] * theta0[2];
+      C[i][j] = cos(thetaNorm) * solidmech::I[i][j] + c1 * theta0Mat_ij + 
+              c2 * theta0[i] * theta0[j];
+    }
+  }
+  storage->getNodePosition(masterNode, masterPos.ptr(), false);
+  for (size_t n = 0; n < slaveNodes.size(); n++) {
+    storage->getNodePosition(slaveNodes[n], slavePos.ptr(), false);
+    Vec<3> pqVec = slavePos - masterPos;
+    // drivatives for C matrix for i row (with respect to j and l)
+    Mat<3,3> cDeriv;
+    for (size_t d = 0; d < dofs.size(); d++) {
+      uint16 i;
+      switch (dofs[d]) {
+        case Dof::UX:
+          i = 0;
+          break;
+        case Dof::UY:
+          i = 1;
+          break;
+        case Dof::UZ:
+          i = 2;
+          break;
+        default:
+          LOG(FATAL) << "which degree of freedom?";
+        }
+        for (uint16 j = 0; j < 3; j++) {
+          for (uint16 l = 0; l < 3; l++) {
+            cDeriv[j][l] = solidmech::I[i][j]*(-c1)*theta0[l] + c3*theta0[l]*(
+                  solidmech::LeviCivita[i][0][j]*theta0[0] +
+                  solidmech::LeviCivita[i][1][j]*theta0[1] +
+                  solidmech::LeviCivita[i][2][j]*theta0[2]) +
+                c1*solidmech::LeviCivita[i][l][j] +
+                c4*theta0[l]*theta0[i]*theta0[j] +
+                c2*(solidmech::I[i][l]*theta0[j] + theta0[i]*solidmech::I[j][l]);
+          }
+        }
+        collection[n*static_cast<uint16> (dofs.size()) + d]->b = w0[i] + (C[i][0] - solidmech::I[i][0])*pqVec[0] +
+              (C[i][1] - solidmech::I[i][1])*pqVec[1] + (C[i][2] - solidmech::I[i][2])*pqVec[2] -
+              storage->getNodeDofSolution(slaveNodes[n], dofs[d]);
+        list<MpcTerm>::iterator token = collection[n*dofs.size() + d]->eq.begin();
+        token++;
+        token++;
+        // masterNode::ROTX coef
+        token->coef = - (cDeriv[0][0]*pqVec[0] + cDeriv[1][0]*pqVec[1] + cDeriv[2][0]*pqVec[2]);
+        token++;
+        // masterNode::ROTY coef
+        token->coef = - (cDeriv[0][1]*pqVec[0] + cDeriv[1][1]*pqVec[1] + cDeriv[2][1]*pqVec[2]);
+        token++;
+        // masterNode::ROTZ coef
+        token->coef = - (cDeriv[0][2]*pqVec[0] + cDeriv[1][2]*pqVec[1] + cDeriv[2][2]*pqVec[2]);
+    }
+  }
+}
+
+
+} // namespace nla3d

nla3d/Mpc.h

+// This file is a part of nla3d project. For information about authors and
+// licensing go to project's repository on github:
+// https://github.com/dmitryikh/nla3d 
+
+#pragma once
+#include <list>
+#include "sys.h"
+#include "Dof.h"
+
+namespace nla3d {
+
+class FEStorage;
+
+
+class MpcTerm {
+public:
+  MpcTerm() : node(0), node_dof(Dof::UNDEFINED), coef(1.0) {
+
+  }
+
+  MpcTerm(int32 n, Dof::dofType dof, double coef = 1.0) : node(n), node_dof(dof), coef(coef) {
+
+  }
+
+  int32 node = 0;
+  Dof::dofType node_dof = Dof::UNDEFINED;
+  double coef = 0.0;
+};
+
+
+class Mpc {
+  public:
+    void print (std::ostream& out);
+    std::list<MpcTerm> eq;
+    double b; // rhs of MPC eq.
+    uint32 eqNum; // number of equation in global assembly
+};
+
+
+class MpcCollection {
+  public:
+    virtual void update() = 0;
+    virtual void pre() = 0;
+    void printEquations(std::ostream& out);
+    void registerMpcsInStorage();
+    std::vector<Mpc*> collection;
+    FEStorage* storage;
+};
+
+
+class RigidBodyMpc : public MpcCollection {
+  public:
+    RigidBodyMpc();
+    ~RigidBodyMpc();
+    uint32 masterNode;
+
+    void pre ();
+    void update ();
+    std::vector<uint32> slaveNodes;
+    std::vector<Dof::dofType> dofs;
+};
+
+
+class fixBC {
+public:
+  fixBC () { }
+  fixBC (int32 n, Dof::dofType dof, double val = 0.0) : node(n), node_dof(dof), value(val)
+  { }
+  int32 node = 0;
+  Dof::dofType node_dof = Dof::UNDEFINED;
+  double value = 0.0;
+};
+
+
+class loadBC {
+public:
+  loadBC () : node(0), node_dof(Dof::UNDEFINED), value(0.0)
+  { }
+  loadBC (int32 n, Dof::dofType dof, double val = 0.0) : node(n), node_dof(dof), value(val)
+  { }
+  int32 node = 0;
+  Dof::dofType node_dof = Dof::UNDEFINED;
+  double value = 0.0;
+};
+
+} // namespace nal3d

nla3d/Node.cpp

+// This file is a part of nla3d project. For information about authors and
+// licensing go to project's repository on github:
+// https://github.com/dmitryikh/nla3d 
+
+#include "Node.h"
+
+namespace nla3d {
+
+
+void Node::display (uint32 nn) {
+  LOG(INFO) << "N " << nn << " : " << pos;
+}
+
+
+std::string Node::toString() {
+  std::string str;
+  char buff[100];
+  sprintf_s(buff,100,"%f %f %f",pos[0], pos[1], pos[2]);
+  str+=buff;
+  return str; //TODO: do it easy
+}
+
+}

nla3d/Node.h

+// This file is a part of nla3d project. For information about authors and
+// licensing go to project's repository on github:
+// https://github.com/dmitryikh/nla3d 
+
+#pragma once
+#include "sys.h"
+#include "math/Vec.h"
+#include "Dof.h"
+
+namespace nla3d {
+
+//class Node represents spatial 3D node
+class Node {
+  public:
+    //in-out operation: (rudiment actually..)
+    void display (uint32 nn);
+    std::string toString();
+
+    math::Vec<3> pos;
+
+    friend class Element;
+};
+
+} // namespace nla3d 

nla3d/PostProcessor.cpp

+// This file is a part of nla3d project. For information about authors and
+// licensing go to project's repository on github:
+// https://github.com/dmitryikh/nla3d 
+
+#include "PostProcessor.h"
+
+namespace nla3d {
+
+PostProcessor::PostProcessor(FEStorage *st) {
+	assert(st);
+	storage = st;
+	active = false;
+	failed = false;
+}
+
+std::string PostProcessor::getStatus () {
+	char buf[300];
+	sprintf_s(buf, 200, "PostProcessor No %d: %s\n\tActive - %s, Failed - %s",nPost_proc, name.c_str(),active?"true":"false",active?"true":"false");
+	return std::string(buf);
+}
+
+} // namespace nla3d 

nla3d/PostProcessor.h

+// This file is a part of nla3d project. For information about authors and
+// licensing go to project's repository on github:
+// https://github.com/dmitryikh/nla3d 
+
+#pragma once
+#include "sys.h"
+#include "FEStorage.h"
+
+namespace nla3d {
+
+class FEStorage;
+class FESolver;
+
+//Data_Processor
+class PostProcessor {
+public:
+	PostProcessor(FEStorage *st);
+	~PostProcessor() { };
+	virtual void pre ()=0;
+	virtual void process (uint16 curLoadstep)=0;
+	virtual void post (uint16 curLoadstep)=0;
+	std::string getStatus ();
+	uint16 getnPost_num () {
+		return nPost_proc;
+	}
+	void setActive (bool act) {
+		if (failed) {
+			LOG(WARNING) << "PostProcessor " << nPost_proc << " (" << name << "): can't set active,"
+          << "because post_proc has been already failed";
+			return;
+		}
+		active = act;
+	}
+	bool getActive () {
+		return active;
+	}
+
+	uint16 getProc_num () {
+		return nPost_proc;
+	}
+
+	friend class FEStorage;
+	friend class FESolver;
+protected:
+	FEStorage *storage;
+	uint16 nPost_proc;
+	std::string name;
+	bool active;
+	bool failed;
+};
+
+} // namespace nla3d

nla3d/ReactionProcessor.cpp

+// This file is a part of nla3d project. For information about authors and
+// licensing go to project's repository on github:
+// https://github.com/dmitryikh/nla3d 
+
+#include "ReactionProcessor.h"
+#include "Node.h"
+
+namespace nla3d {
+
+ReactionProcessor::ReactionProcessor(FEStorage *st) : PostProcessor(st) {
+	name ="ReactionProcessor";
+}
+
+ReactionProcessor::ReactionProcessor(FEStorage *st, std::string _filename) : PostProcessor(st) {
+	name ="ReactionProcessor";
+  filename = _filename;
+}
+
+void ReactionProcessor::pre() {
+	if (nodes.size() == 0) {
+    LOG(WARNING) << "Can't work. No nodes. Processor name = " << name;
+    return;
+  }
+
+  // define Dofs by the following rule:
+  // 1. If dofs is not empty than user alreade define which dof use to calculate reaction
+  // 2. If dofs.size() == 0 , select dof thich are constrained in current node set.
+  // if different nodes has different constrained dofs - error
+  if (dofs.size() == 0) {
+    for (uint16 d = 0; d < Dof::numberOfDofTypes; d++) {
+      Dof::dofType t = static_cast<Dof::dofType> (d);
+      if (storage->isNodeDofUsed(nodes[0], t) && storage->getNodeDof(nodes[0], t)->isConstrained) {
+        dofs.push_back(t);
+      } 
+    }
+    if (dofs.size() == 0) {
+      LOG(WARNING) << "Can't select dofs for reactions (they are unconstrained)";
+    } 
+    // check that others nodes have the same constrained dofs
+    bool sameDofsConstrained = true;
+    uint32 n = 0;
+    while (sameDofsConstrained && n < nodes.size()) {
+      for (uint16 d = 0; d < dofs.size(); d++) {
+        if (!storage->isNodeDofUsed(nodes[n], dofs[d]) && !storage->getNodeDof(nodes[n], dofs[d])->isConstrained) {
+          LOG(FATAL) << "Different dofs are constrained in the node set. Autochoosing of dofs is failed";
+          sameDofsConstrained = false;
+          break;
+        }
+      }
+      n++;
+    }
+  } else {
+    // in case of user defined dofs just check that this dofs are used in the FE
+    bool dofsUsed = true;
+    uint32 n = 0;
+    while (dofsUsed && n < nodes.size()) {
+      for (uint16 d = 0; d < dofs.size(); d++) {
+        if (!storage->isNodeDofUsed(nodes[n], dofs[d])) {
+          LOG(ERROR) << "some dofs are not used in FE calculations";
+          dofsUsed = false;
+          break;
+        }
+      }
+      n++;
+    }
+  }
+  if (filename.length() > 0) {
+    std::ofstream file(filename.c_str(),std::ios::trunc);
+    if (!file) {
+      LOG(WARNING) << "Can't create a file with name " << filename;
+      return;
+    }
+    for (uint16 d = 0; d < dofs.size(); d++) {
+      file << "\t" << Dof::dofTypeLabels[dofs[d]];
+    }
+    file << std::endl;
+    file << "0";
+    for (uint16 d = 0; d < dofs.size(); d++) {
+      file << "\t" << 0.0;
+    }
+    file << std::endl;
+    file.close();
+  }
+  sumOfDofsReactions.assign(dofs.size(), std::vector<double> ());
+  for (uint16 d = 0; d < dofs.size(); d++) {
+    //sumOfDofsReactions[d].reserve(qLoadstep+1);
+    sumOfDofsReactions[d].push_back(0.0);
+  }
+}
+
+void ReactionProcessor::process (uint16 curLoadstep) {
+  std::vector<double> reactions;
+  reactions.assign(dofs.size(), 0.0);
+
+	for (uint32 n = 0; n < nodes.size(); n++) {
+    for (uint16 d = 0; d < dofs.size(); d++) {
+      reactions[d] += storage->getReaction(nodes[n],dofs[d]);
+    }
+  }
+
+  if (filename.length() > 0) {
+    std::ofstream file(filename.c_str(),std::ios::app);
+    if (!file) {
+      LOG(WARNING) << "Can't open a file with name " << filename;
+      return;
+    }
+    file << curLoadstep;
+    for (uint16 d = 0; d < dofs.size(); d++) {
+      file << "\t" << reactions[d];
+    }
+    file << std::endl;
+    file.close();
+  }
+  for (uint16 d = 0; d < dofs.size(); d++) {
+    sumOfDofsReactions[d].push_back(reactions[d]);
+  }
+}
+
+
+void ReactionProcessor::post (uint16 curLoadstep) {
+}
+
+
+std::vector<double> ReactionProcessor::getReactions (Dof::dofType dof) {
+  for (uint16 d = 0; d < dofs.size(); d++) {
+    if (dofs[d] == dof) {
+      return sumOfDofsReactions[d];
+    }
+  }
+  LOG(ERROR) << "There is not results for DoF = " << Dof::dofTypeLabels[dof];
+  return std::vector<double> ();
+}
+
+} // namespace nla3d

nla3d/ReactionProcessor.h

+// This file is a part of nla3d project. For information about authors and
+// licensing go to project's repository on github:
+// https://github.com/dmitryikh/nla3d 
+
+#pragma once
+#include "sys.h"
+#include "PostProcessor.h"
+#include "FEStorage.h"
+
+namespace nla3d {
+
+class ReactionProcessor : public PostProcessor {
+public:
+	ReactionProcessor(FEStorage *st);
+	ReactionProcessor(FEStorage *st, std::string _filename);
+	virtual ~ReactionProcessor() { };
+
+	virtual void pre ();
+	virtual void process (uint16 curLoadstep);
+	virtual void post (uint16 curLoadstep);
+
+  std::vector<double> getReactions (Dof::dofType dof);
+	
+	std::vector<uint32> nodes;
+	std::vector<Dof::dofType> dofs;
+protected:
+  std::string filename;
+	std::vector<std::vector<double> > sumOfDofsReactions;
+};
+
+} // namespace nla3d

nla3d/VtkProcessor.h

+// This file is a part of nla3d project. For information about authors and
+// licensing go to project's repository on github:
+// https://github.com/dmitryikh/nla3d 
+
+#pragma once
+#include <set>
+#include <map>
+#include "sys.h"
+#include "PostProcessor.h"
+#include "FEStorage.h"
+#include "elements/query.h"
+
+
+namespace nla3d {
+
+// VtkProcessor is dedicated to dump mesh and results into *.vtk legacy textual format. On pre()
+// callback filename0.vtk is written which contains initial state of FE problem before solution. On
+// every solution step process(...) callback is called which writes filename%d.vtk with solution
+// data relevant to this step.  Later this banch of files could be used to visualize solution fields
+// in external program (ex.  Paraview).
+//
+// By default, VtkProcessor write DoF fields (displacements, reactions for structural problems and
+// temperatures for thermal). But by registering element results by registerResult(...) one can add
+// addition fields (scalar, vector or tensor ones) to be written into vtk files.
+
+class VtkProcessor : public PostProcessor {
+public:
+	VtkProcessor(FEStorage *st, std::string _fileName);
+	virtual ~VtkProcessor();
+	virtual void pre();
+	virtual void process(uint16 curLoadstep);
+	virtual void post(uint16 curLoadstep);
+
+  // One can use writeAllResults(true) in order to ask VtkProcess to determine which query codes is
+  // relevant to FEStorage's elements automatically. As results, all accessible element results will
+  // be written into vtk files.
+  void writeAllResults(bool write = true);
+
+  // The methods below can be used to manually set query codes for elements. Every query will be add
+  // addition cell data fields into vtk file.
+  bool registerResults(std::initializer_list<scalarQuery> queries);
+  bool registerResults(std::initializer_list<vectorQuery> queries);
+  bool registerResults(std::initializer_list<tensorQuery> queries);
+
+private:
+	std::string file_name;
+
+	void write_header(std::ofstream &file);
+	void write_header_binary();
+  // write geometry(mesh) into the vtk's `file`. If `def == true` then write node coordinates in
+  // deformed state
+	void write_geometry(std::ofstream &file, bool def = false);
+	void write_geometry_binary(bool def = false);
+	void write_point_data(std::ofstream &file);
+	void write_point_data_binary();
+	void write_cell_data(std::ofstream &file);
+	void write_cell_data_binary();
+
+  void writeScalar(std::ofstream &file, const char* name, std::vector<double>& data);
+  void writeScalar_binary(const char* name, std::vector<double>& data);
+  void writeTensor(std::ofstream &file, const char* name, std::vector<math::MatSym<3> >& data);
+  void writeTensor_binary(const char* name, std::vector<math::MatSym<3> >& data);
+  void writeVector(std::ofstream &file, const char* name, std::vector<math::Vec<3> >& data);
+  void writeVector_binary(const char* name, std::vector<math::Vec<3> >& data);
+  void write_vecres(std::ofstream &file);
+
+  void revealAllResults();
+
+  bool useDisplacementsDofs = false;
+  bool _writeAllResults = false;
+
+  std::set<Dof::dofType> nodalDofTypes;
+  std::set<Dof::dofType> elementDofTypes;
+
+  std::set<scalarQuery> cellScalarQueries;
+  std::set<vectorQuery> cellVectorQueries;
+  std::set<tensorQuery> cellTensorQueries;
+};
+
+} // namespace nla3d

nla3d/elements/ElementFactory.cpp

+// This file is a part of nla3d project. For information about authors and
+// licensing go to project's repository on github:
+// https://github.com/dmitryikh/nla3d 
+
+#include "elements/ElementFactory.h"
+#include "elements/PLANE41.h"
+#include "elements/SOLID81.h"
+#include "elements/TRUSS3.h"
+#include "elements/TETRA0.h"
+#include "elements/TETRA1.h"
+#include "elements/QUADTH.h"
+#include "elements/TRIANGLE4.h"
+
+namespace nla3d {
+
+
+ElementType ElementFactory::elName2elType (std::string elName) {
+  for (uint16 i = 0; i < (uint16) ElementType::UNDEFINED; i++) {
+    if (elName.compare(elTypeLabels[i]) == 0) {
+      return (ElementType) i;
+    }
+  }
+  return ElementType::UNDEFINED;
+}
+
+
+void ElementFactory::createElements (ElementType elId, const uint32 n, std::vector<Element*>& ptr) {
+  if (elId == ElementType::UNDEFINED) {
+    LOG(FATAL) << "Element type is undefined";
+  }
+  if (n == 0) {
+    return;
+  }
+  // reserve memory in the vector for newly created elements
+  ptr.reserve(ptr.size() + n);
+
+  switch (elId) {
+      case ElementType::PLANE41:
+        for (uint32 i = 0; i < n; i++) {
+          ptr.push_back(new ElementPLANE41());
+        }
+        break;
+      case ElementType::SOLID81:
+        for (uint32 i = 0; i < n; i++) {
+          ptr.push_back(new ElementSOLID81());
+        }
+        break;
+      case ElementType::TRUSS3:
+        for (uint32 i = 0; i < n; i++) {
+          ptr.push_back(new ElementTRUSS3());
+        }
+        break;
+      case ElementType::TRIANGLE4:
+        for (uint32 i = 0; i < n; i++) {
+          ptr.push_back(new ElementTRIANGLE4());
+        }
+        break;
+      case ElementType::TETRA0:
+        for (uint32 i = 0; i < n; i++) {
+          ptr.push_back(new ElementTETRA0());
+        }
+        break;
+      case ElementType::TETRA1:
+        for (uint32 i = 0; i < n; i++) {
+          ptr.push_back(new ElementTETRA1());
+        }
+        break;
+      case ElementType::QUADTH:
+        for (uint32 i = 0; i < n; i++) {
+          ptr.push_back(new ElementQUADTH());
+        }
+        break;
+      case ElementType::SurfaceLINETH:
+        for (uint32 i = 0; i < n; i++) {
+          ptr.push_back(new SurfaceLINETH());
+        }
+        break;
+      default:
+        LOG(ERROR) << "Don't have an element with id " << (uint16) elId;
+    }
+}
+
+} // namespace nla3d

nla3d/elements/ElementFactory.h

+// This file is a part of nla3d project. For information about authors and
+// licensing go to project's repository on github:
+// https://github.com/dmitryikh/nla3d 
+
+#pragma once
+#include "sys.h"
+
+namespace nla3d {
+class Element;
+
+class ElementFactory {
+  public:
+
+    static ElementType elName2elType (std::string elName); 
+    static void createElements (ElementType elId, const uint32 n, std::vector<Element*>& ptr); 
+};
+
+} // namespace nla3d

nla3d/elements/PLANE41.h

+// This file is a part of nla3d project. For information about authors and
+// licensing go to project's repository on github:
+// https://github.com/dmitryikh/nla3d 
+
+#pragma once
+#include "elements/element.h"
+#include "elements/isoparametric.h"
+#include "FEStorage.h"
+#include "solidmech.h"
+
+namespace nla3d {
+
+//4-node 2D QUAD nonlinear element based on mixed approach
+class ElementPLANE41 : public ElementIsoParamQUAD {
+  public:
+    ElementPLANE41 () {
+      intOrder = 2;
+      type = ElementType::PLANE41;
+    }
+
+    ElementPLANE41 (const ElementPLANE41& from) {
+      operator=(from);
+    }
+
+    //solving procedures
+    void pre();
+    void buildK();
+    void update();
+    math::Mat<3,8> make_B (uint16 nPoint);  //функция создает линейную матрицу [B]
+    math::Mat<4,8> make_Bomega (uint16 nPoint); //функция создает линейную матрицу [Bomega]
+
+    //postproc procedures
+    bool getScalar(double* scalar, scalarQuery code, uint16 gp, const double scale);
+    bool getVector(math::Vec<3>* vector, vectorQuery code, uint16 gp, const double scale);
+    bool getTensor(math::MatSym<3>* tensor, tensorQuery code, uint16 gp, const double scale);
+
+    // internal element data
+    // S[0] - Sx  S[1] - Sy S[2] - Sxy
+    std::vector<math::Vec<3> > S; //S[номер т. интегр.][номер напряжения] - напряжения Пиолы-Кирхгоффа
+    // C[0] - C11 C[1] - C22  C[2] - C12
+    std::vector<math::Vec<3> > C; //C[номер т. интегр.][номер деформ.] - компоненты матрицы меры деформации
+    // O[0] - dU/dx O[1] - dU/dy  O[2] - dV/dx  O[3] - dV/dy
+    std::vector<math::Vec<4> > O; //S[номер т. интегр.][номер омеги]
+
+    // addition data
+    static const solidmech::tensorComponents components[3];
+    static const uint16 num_components;
+
+
+    template <uint16 el_dofs_num>
+    void assembleK(const math::Mat<el_dofs_num,el_dofs_num> &Ke, const math::Vec<el_dofs_num> &Qe);
+};
+
+template <uint16 el_dofs_num>
+void ElementPLANE41::assembleK(const math::Mat<el_dofs_num,el_dofs_num> &Ke, const math::Vec<el_dofs_num> &Qe)
+{
+  uint16 dim = 2;
+  uint16 eds = 8;
+  uint16 eldofs = 1;
+  Dof::dofType nodeDofVec[] = {Dof::UX, Dof::UY};
+  Dof::dofType elementDofVec[] = {Dof::HYDRO_PRESSURE};
+  assert(getNNodes() * dim + eldofs == el_dofs_num);
+  // upper diagonal process for nodal dofs
+  for (uint16 i=0; i < getNNodes(); i++)
+    for (uint16 j=i; j < getNNodes(); j++)
+      for (uint16 di=0; di < dim; di++)
+        for (uint16 dj=0; dj < dim; dj++)
+          if ((i==j) && (dj>di)) continue;
+          else
+            storage->addValueK(nodes[i], nodeDofVec[di],nodes[j],nodeDofVec[dj], Ke[i*dim+di][j*dim+dj]);
+
+  //upper diagonal process for nodes-el dofs
+  for (uint16 i=0; i < getNNodes(); i++)
+    for(uint16 di=0; di < dim; di++) {
+      uint32 noEq = storage->getNodeDofEqNumber(nodes[i], nodeDofVec[di]);
+      for (uint16 dj=0; dj < eldofs; dj++) {
+        uint32 elEq = storage->getElementDofEqNumber(getElNum(), elementDofVec[dj]);
+        storage->addValueK(noEq, elEq, Ke[i*dim+di][eds+dj]);
+      }
+    }
+  //upper diagonal process for el-el dofs
+  for (uint16 di=0; di < eldofs; di++)
+    for (uint16 dj=di; dj < eldofs; dj++) {
+      uint32 elEqi = storage->getElementDofEqNumber(getElNum(), elementDofVec[di]);
+      uint32 elEqj = storage->getElementDofEqNumber(getElNum(), elementDofVec[dj]);
+      storage->addValueK(elEqi, elEqj, Ke[eds+di][eds+dj]);
+    }
+
+  for (uint16 i=0; i < getNNodes(); i++)
+    for (uint16 di=0; di < dim; di++)
+      storage->addValueF(nodes[i], nodeDofVec[di], Qe[i*dim+di]);
+
+  for (uint16 di=0; di < eldofs; di++) {
+    uint32 elEq = storage->getElementDofEqNumber(getElNum(), elementDofVec[di]);
+    storage->addValueF(elEq, Qe[eds+di]);
+  }
+}
+
+} // namespace nla3d 

nla3d/elements/QUADTH.cpp

+// This file is a part of nla3d project. For information about authors and
+// licensing go to project's repository on github:
+// https://github.com/dmitryikh/nla3d 
+
+#include "elements/QUADTH.h"
+
+namespace nla3d {
+using namespace math;
+
+
+void ElementQUADTH::pre() {
+  if (det.size()==0) {
+    makeJacob();
+  }
+
+  // register element equations
+  for (uint16 i = 0; i < getNNodes(); i++) {
+    storage->addNodeDof(getNodeNumber(i), {Dof::TEMP});
+  }
+}
+
+void ElementQUADTH::buildK() {
+  MatSym<4> Ke;  // element stiff. matrix
+  Ke.zero();
+
+  Vec<4> Fe; // rhs of element equations
+  Fe.zero();
+
+  MatSym<2> mat_k;
+  mat_k.zero();
+  mat_k.comp(0,0) = k;
+  mat_k.comp(1,1) = k;
+
+  // build Ke
+  double dWt; //Gaussian quadrature weight
+  for (uint16 np=0; np < nOfIntPoints(); np++) {
+    dWt = intWeight(np);
+    Mat<2,4> matB = make_B(np);
+    matBTDBprod(matB, mat_k, dWt, Ke);
+
+    // for volume flux
+    if (volFlux != 0.0) {
+      Fe += formFunc(np) * (volFlux * dWt);
+    }
+  }// loop over integration points
+
+  assembleK(Ke, Fe, {Dof::TEMP});
+}
+
+
+void ElementQUADTH::buildC() {
+  MatSym<4> Ce;  // element stiff. matrix
+  Ce.zero();
+
+
+  // build Ke
+  double dWt; //Gaussian quadrature weight
+  for (uint16 np=0; np < nOfIntPoints(); np++) {
+    dWt = intWeight(np);
+    Vec<4> ff = formFunc(np);
+    for (uint16 i = 0; i < 4; i++)
+      for (uint16 j = i; j < 4; j++) 
+        Ce.comp(i, j) += ff[i] * ff[j] * rho * c * dWt;
+  }// loop over integration points
+
+  assembleC(Ce, {Dof::TEMP});
+}
+
+inline Mat<2,4> ElementQUADTH::make_B(uint16 np) {
+  return Mat<2,4>(NiXj[np][0][0], NiXj[np][1][0], NiXj[np][2][0], NiXj[np][3][0],
+                  NiXj[np][0][1], NiXj[np][1][1], NiXj[np][2][1], NiXj[np][3][1]);
+}
+
+
+void ElementQUADTH::update() {
+
+}
+
+void SurfaceLINETH::pre() {
+  makeJacob();
+
+  // register element equations
+  for (uint16 i = 0; i < getNNodes(); i++) {
+    storage->addNodeDof(getNodeNumber(i), {Dof::TEMP});
+  }
+
+  // if environment temperature isn't defined for node 2, we suppose that they are equal
+  if (etemp[1] == 0.0) {
+    etemp[1] = etemp[0];
+  }
+}
+
+
+void SurfaceLINETH::buildK() {
+  MatSym<2> Ke;  // element stiff. matrix
+  Ke.zero();
+
+  Vec<2> Fe; // rhs of element equations
+  Fe.zero();
+
+  double dWt; //Gaussian quadrature weight
+
+  // flux over surface (line)
+  if (flux != 0.0) {
+    for (uint16 np = 0; np < nOfIntPoints(); np++) {
+      dWt = intWeight(np);
+      Vec<2> Ns = formFunc(np);
+      Fe += Ns * (flux * dWt);
+    }
+  }
+
+  // convection over surface (line)
+  if (htc != 0.0) {
+    for (uint16 np = 0; np < nOfIntPoints(); np++) {
+      dWt = intWeight(np);
+      Vec<2> Ns = formFunc(np);
+      MatSym<2> tmp;
+      tmp.comp(0,0) = Ns[0] * Ns[0] * htc * dWt;
+      tmp.comp(0,1) = Ns[0] * Ns[1] * htc * dWt;
+      tmp.comp(1,0) = tmp.comp(0,1);
+      tmp.comp(1,1) = Ns[1] * Ns[1] * htc * dWt;
+      //matAATprod(Ns, htc * dWt, tmp);
+      Ke += tmp;
+      matBVprod(tmp, etemp, 1.0, Fe);
+    }
+  }
+
+  assembleK(Ke, Fe, {Dof::TEMP});
+}
+
+void SurfaceLINETH::update() {
+
+}
+
+
+} // namespace nla3d

nla3d/elements/QUADTH.h

+// This file is a part of nla3d project. For information about authors and
+// licensing go to project's repository on github:
+// https://github.com/dmitryikh/nla3d 
+
+#pragma once
+#include "elements/element.h"
+#include "elements/isoparametric.h"
+#include "FEStorage.h"
+#include "solidmech.h"
+
+namespace nla3d {
+
+//4-node 2D QUAD  element for steady thermal analysis
+class ElementQUADTH : public ElementIsoParamQUAD {
+  public:
+    ElementQUADTH () {
+      intOrder = 2;
+      type = ElementType::QUADTH;
+    }
+
+    //solving procedures
+    void pre();
+    void buildK();
+    void buildC();
+    void buildM() { };
+    void update();
+    math::Mat<2,4> make_B (uint16 nPoint);  // make derivatives matrix
+
+    // conductivity coef ( W/(K m), for example)
+    double k = 0.0;
+    double rho = 0.0; // density
+    double c = 0.0; // capacity
+
+    // volume flux
+    double volFlux = 0.0;
+};
+
+class SurfaceLINETH : public ElementIsoParamLINE {
+  public:
+    SurfaceLINETH () {
+      intOrder = 2;
+      type = ElementType::SurfaceLINETH;
+    }
+
+    //solving procedures
+    void pre();
+    void buildK();
+    void buildC() { };
+    void buildM() { };
+    void update();
+
+    double flux = 0.0;
+    double htc = 0.0;
+    math::Vec<2> etemp = {0.0, 0.0};
+};
+} // nla3d namespace