diff --git a/src/FEReaders.cpp b/src/FEReaders.cpp index 05bd890..429455b 100644 --- a/src/FEReaders.cpp +++ b/src/FEReaders.cpp @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #include "FEReaders.h" #include @@ -27,14 +27,14 @@ void MeshData::clear() { nodesPos.clear(); loadBcs.clear(); fixBcs.clear(); - mpcs.clear(); + mpcs.clear(); feComps.clear(); } void MeshData::compressNumbers() { // numbers should start from 1 - if (nodesNumbers.size() == *std::max_element(nodesNumbers.begin(), nodesNumbers.end())) { + if (nodesNumbers.size() == *std::max_element(nodesNumbers.begin(), nodesNumbers.end())) { // nothing to do return; } @@ -77,7 +77,7 @@ void MeshData::compressNumbers() { if (comp.type == FEComponent::NODES) { for (auto n : comp.list) { n = old2new[n]; - } + } } } } @@ -163,7 +163,7 @@ std::vector ssplit(const std::string& line, const std::vector& return vv; } - + bool iequals(const string& a, const string& b) { unsigned int sz = a.size(); if (b.size() != sz) @@ -175,7 +175,7 @@ bool iequals(const string& a, const string& b) { } // taken from -// http://stackoverflow.com/questions/6089231/getting-std-ifstream-to-handle-lf-cr-and-crlf/6089413#6089413 +// http://stackoverflow.com/questions/6089231/getting-std-ifstream-to-handle-lf-cr-and-crlf/6089413#6089413 std::istream& getLine(std::istream& is, std::string& t) { t.clear(); @@ -273,7 +273,7 @@ bool readCdbFile(std::string filename, MeshData& md) { // current element type number. We keep the current element type number while proceed the apld // file. Currently, this info is not used, but in the past it was used to determine element type - // while parsing EBLOCK blocks. + // while parsing EBLOCK blocks. uint32 type = 0; Tokenizer t; @@ -301,7 +301,7 @@ bool readCdbFile(std::string filename, MeshData& md) { // 1 2.000000000E+000 6.000000000E+000 0.000000000E+000 // 2 0.000000000E+000 6.000000000E+000 0.000000000E+000 // - // from Ansys WB APDL 17.1: + // from Ansys WB APDL 17.1: //NBLOCK,6,SOLID, 341, 341 //(3i9,6e21.13e3) // 1 0 0 0.0000000000000E+000 4.0000000000000E+000 @@ -352,7 +352,7 @@ bool readCdbFile(std::string filename, MeshData& md) { md.nodesPos.push_back(pos); } }//NBLOCK - else if (iequals(t.tokens[0], "EBLOCK")) { + else if (iequals(t.tokens[0], "EBLOCK")) { // Example of record: //EBLOCK,19,SOLID, 7024, 7024 //(19i9) @@ -445,7 +445,7 @@ bool readCdbFile(std::string filename, MeshData& md) { st = 5; } vector enodes; - + if (v.size() != st + nNodes) { LOG(FATAL) << "Not enought fields to read element nodes"; } @@ -484,8 +484,8 @@ bool readCdbFile(std::string filename, MeshData& md) { else if (iequals(t.tokens[0], "CE")) { // CE - is a command to declare MPC equation // How MPC looks like this in cdb file: - //CE,R5.0,DEFI, 2, 1, 0.00000000 - //CE,R5.0,NODE, 1700,UX , 1.00000000 , 1700,UZ , 1.00000000 + //CE,R5.0,DEFI, 2, 1, 0.00000000 + //CE,R5.0,NODE, 1700,UX , 1.00000000 , 1700,UZ , 1.00000000 // TODO: Mpc is dynamicaly allocated, but MeshData won't free this memory (this should be done // in FEStorage). This is potential memory leak Mpc* mpc = new Mpc(); diff --git a/src/FEReaders.h b/src/FEReaders.h index 0e07f3c..e21be32 100644 --- a/src/FEReaders.h +++ b/src/FEReaders.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once #include @@ -38,7 +38,7 @@ class MeshData { std::vector loadBcs; std::vector fixBcs; - std::vector mpcs; + std::vector mpcs; std::map feComps; diff --git a/src/lib/Dof.cpp b/src/lib/Dof.cpp index 95479ba..ba8b8a1 100644 --- a/src/lib/Dof.cpp +++ b/src/lib/Dof.cpp @@ -1,3 +1,7 @@ +// 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 "Dof.h" namespace nla3d { diff --git a/src/lib/Dof.h b/src/lib/Dof.h index d1dc545..a772326 100644 --- a/src/lib/Dof.h +++ b/src/lib/Dof.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once #include "sys.h" @@ -8,7 +8,7 @@ namespace nla3d { class DofCollection; -// class for Degree of Freedom informations +// class for Degree of Freedom informations // (is it fixed boundary condition, what its number in equation system, ..) class Dof { public: @@ -65,10 +65,10 @@ class DofCollection { uint32 numberOfUsedDofs = 0; uint32 numberOfEntities = 0; - // vector of Dof objects + // vector of Dof objects std::vector 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 + // Dof for entity n will be located from dofPos[n-1] included to dofPos[n] excluded std::vector dofPos; // set of unique dofs used in collection std::set uniqueDofTypes; @@ -95,4 +95,4 @@ inline std::pair::iterator, } -} // namespace nla3d +} // namespace nla3d diff --git a/src/lib/FEComponent.cpp b/src/lib/FEComponent.cpp index 49adba2..87a4c9a 100644 --- a/src/lib/FEComponent.cpp +++ b/src/lib/FEComponent.cpp @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #include "FEComponent.h" @@ -14,7 +14,7 @@ FEComponent::FEComponent () { } FEComponent::~FEComponent () { list.clear(); -} +} FEComponent::typeOfComponent FEComponent::typeFromString(const std::string& typeName) { for (size_t i = 1; i < LAST; i++) { diff --git a/src/lib/FEComponent.h b/src/lib/FEComponent.h index 3f2b4cc..da11318 100644 --- a/src/lib/FEComponent.h +++ b/src/lib/FEComponent.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once #include "sys.h" @@ -32,6 +32,6 @@ class FEComponent { inline MAKE_LOGGABLE(FEComponent, obj, os) { os << obj.name << ": " << obj.list.size() << " " << obj.labelsOfComponent[obj.type]; return os; -} +} -} // namespace nla3d +} // namespace nla3d diff --git a/src/lib/FESolver.cpp b/src/lib/FESolver.cpp index 148e005..95b9429 100644 --- a/src/lib/FESolver.cpp +++ b/src/lib/FESolver.cpp @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #include "FESolver.h" @@ -27,7 +27,7 @@ uint16 TimeControl::getCurrentEquilibriumStep() { uint16 TimeControl::getTotalNumberOfEquilibriumSteps() { - return totalNumberOfEquilibriumSteps; + return totalNumberOfEquilibriumSteps; } @@ -37,7 +37,7 @@ bool TimeControl::nextStep(double delta) { convergedTimeInstances.push_back(currentTime); } if (currentTime >= endTime) { - return false; + return false; } // if delta bigger than endTime-currentTime if (currentTime + delta >= endTime) { @@ -58,7 +58,7 @@ bool TimeControl::nextStep(double delta) { void TimeControl::nextEquilibriumStep() { currentEquilibriumStep++; - totalNumberOfEquilibriumSteps++; + totalNumberOfEquilibriumSteps++; LOG(INFO) << "***** Equilibrium iteration = " << currentEquilibriumStep << ", Cumulutive iterations = " << totalNumberOfEquilibriumSteps; } @@ -229,7 +229,7 @@ void FESolver::setConstrainedDofs() { void FESolver::applyBoundaryConditions(double time) { TIMED_SCOPE(t, "applyBoundaryConditions"); LOG(INFO) << "Applying boundary conditions.. (" << loads.size() << " nodal loads and " - << fixs.size() << "nodal fixations)"; + << fixs.size() << "nodal fixations)"; // fill nodal loads for (auto load : loads) { @@ -340,7 +340,7 @@ void LinearFESolver::solve () { CHECK_NOTNULL(storage); CHECK_NOTNULL(eqSolver); - // setup matrix properties for EquationSolver + // setup matrix properties for EquationSolver eqSolver->setSymmetric(true); eqSolver->setPositive(false); @@ -410,7 +410,7 @@ void NonlinearFESolver::solve() { CHECK_NOTNULL(storage); CHECK_NOTNULL(eqSolver); - // setup matrix properties for EquationSolver + // setup matrix properties for EquationSolver eqSolver->setSymmetric(true); eqSolver->setPositive(false); @@ -474,16 +474,16 @@ void NonlinearFESolver::solve() { currentCriteria = calculateCriteria(deltaUs); // TODO: 1. It seems that currentCriteria is already normalized in calculateCriteria(). we - // need to compare currentCriteria with 1.0 + // need to compare currentCriteria with 1.0 // TODO: 2. Current convergence criteria is not good. We need to introduce equilibrium balance - // criteria too along with kinematic one. + // criteria too along with kinematic one. if (currentCriteria < convergenceCriteria) { converged = true; break; } LOG_IF(currentCriteria > 1.0e6 || std::isnan(currentCriteria), FATAL) << "The solution is diverged!"; - if (timeControl.getCurrentEquilibriumStep() >= numberOfIterations) + if (timeControl.getCurrentEquilibriumStep() >= numberOfIterations) break; }//iterations @@ -535,7 +535,7 @@ void LinearTransientFESolver::solve() { CHECK_NOTNULL(storage); CHECK_NOTNULL(eqSolver); - // setup matrix properties for EquationSolver + // setup matrix properties for EquationSolver eqSolver->setSymmetric(true); eqSolver->setPositive(false); @@ -558,7 +558,7 @@ void LinearTransientFESolver::solve() { setConstrainedDofs(); storage->assignEquationNumbers(); initSolutionData(); - + vecR.zero(); @@ -582,7 +582,7 @@ void LinearTransientFESolver::solve() { vecDU.zero(); vecDDU.zero(); - + for (size_t i = 0; i < getNumberOfPostProcessors(); i++) { postProcessors[i]->pre(); } diff --git a/src/lib/FESolver.h b/src/lib/FESolver.h index 2f72603..ed9ce8a 100644 --- a/src/lib/FESolver.h +++ b/src/lib/FESolver.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once #include "sys.h" @@ -118,11 +118,11 @@ class FESolver { // this work based on `fixs` array. void setConstrainedDofs(); - // add DoF fixation (constraint) boundary condition + // add DoF fixation (constraint) boundary condition void addFix(int32 n, Dof::dofType dof, const double value = 0.0); - // add DoF load (force) boundary condition + // 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); @@ -141,7 +141,7 @@ class FESolver { std::vector 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. + // references make it easy to operate with important FE data structures. BlockSparseSymMatrix<2>* matK = nullptr; BlockSparseSymMatrix<2>* matC = nullptr; BlockSparseSymMatrix<2>* matM = nullptr; diff --git a/src/lib/FEStorage.cpp b/src/lib/FEStorage.cpp index 211a983..021682f 100644 --- a/src/lib/FEStorage.cpp +++ b/src/lib/FEStorage.cpp @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #include "FEStorage.h" #include "elements/element.h" @@ -14,7 +14,7 @@ FEStorage::FEStorage() { FEStorage::~FEStorage () { - if (material) { + if (material) { delete material; } @@ -24,7 +24,7 @@ FEStorage::~FEStorage () { void FEStorage::assembleGlobalEqMatrices() { - assert(matK); + assert(matK); assert(matK->isCompressed()); TIMED_SCOPE(t, "assembleGlobalEqMatrix"); @@ -120,26 +120,26 @@ double FEStorage::getReaction(uint32 eq) { } Material* FEStorage::getMaterial() { - assert(material); - return material; + assert(material); + return material; } void FEStorage::getElementNodes(uint32 el, Node** node_ptr) { - assert(el <= nElements()); + assert(el <= nElements()); Element* elp = elements[el-1]; - for (uint16 i=0; igetNNodes(); i++) - node_ptr[i] = nodes[elp->getNodeNumber(i)-1]; + for (uint16 i=0; igetNNodes(); i++) + node_ptr[i] = nodes[elp->getNodeNumber(i)-1]; } // prt array always has 3 elements void FEStorage::getNodePosition(uint32 n, double* ptr, bool deformed) { - assert(n > 0 && n <= nNodes()); - for (uint16 i = 0; i < 3; i++) { - ptr[i] = nodes[n-1]->pos[i]; + assert(n > 0 && n <= nNodes()); + for (uint16 i = 0; i < 3; i++) { + ptr[i] = nodes[n-1]->pos[i]; } - if (deformed) { + if (deformed) { if (isNodeDofUsed(n, Dof::UX)) { ptr[0] += getNodeDofSolution(n, Dof::UX); } @@ -151,7 +151,7 @@ void FEStorage::getNodePosition(uint32 n, double* ptr, bool deformed) { if (isNodeDofUsed(n, Dof::UZ)) { ptr[2] += getNodeDofSolution(n, Dof::UZ); } - } + } } @@ -175,7 +175,7 @@ FEComponent* FEStorage::getFEComponent(const std::string& name) { void FEStorage::addMpc (Mpc* mpc) { assert (mpc); assert (mpc->eq.size() > 0); - mpcs.push_back(mpc); + mpcs.push_back(mpc); } @@ -185,15 +185,15 @@ void FEStorage::addMpcCollection (MpcCollection* mpcCol) { } void FEStorage::addFEComponent (FEComponent* comp) { - assert(comp); - feComponents.push_back(comp); + assert(comp); + feComponents.push_back(comp); } void FEStorage::addNode (Node* node) { //Node() fires Vec<3> constructor, thus Node coordinates are (0,0,0) by default //TODO: try-catch of memory overflow - nodes.push_back(node); + nodes.push_back(node); } @@ -222,7 +222,7 @@ void FEStorage::addElement (Element* el) { el->storage = this; el->elNum = nElements() + 1; //TODO: try-catch of memory overflow - elements.push_back(el); + elements.push_back(el); } @@ -231,7 +231,7 @@ std::vector FEStorage::createElements(uint32 _en, ElementType elType) { std::vector newIndexes; newIndexes.reserve(_en); uint32 nextNumber = elements.size() + 1; - ElementFactory::createElements(elType, _en, elements); + ElementFactory::createElements(elType, _en, elements); for (uint32 i = nextNumber; i <= elements.size(); i++) { //access elNum protected values as friend elements[i - 1]->elNum = i; @@ -243,8 +243,8 @@ std::vector FEStorage::createElements(uint32 _en, ElementType elType) { void FEStorage::deleteMesh() { - deleteElements(); - deleteNodes(); + deleteElements(); + deleteNodes(); deleteMpcs(); deleteMpcCollections(); deleteFeComponents(); @@ -301,9 +301,9 @@ void FEStorage::deleteDofArrays() { void FEStorage::deleteSolutionData() { - _nDofs = 0; + _nDofs = 0; _nConstrainedDofs = 0; - _nUnknownDofs = 0; + _nUnknownDofs = 0; if (matK) { delete matK; @@ -317,7 +317,7 @@ void FEStorage::deleteSolutionData() { delete matM; matM = nullptr; } - + deleteDofArrays(); vecU.clear(); @@ -337,14 +337,14 @@ void FEStorage::listFEComponents() { void FEStorage::initDofs() { TIMED_SCOPE(t, "initDofs"); - if (nNodes() < 1 && nElements() < 1) { - LOG(FATAL) << "No any nodes or elements"; - } - + if (nNodes() < 1 && nElements() < 1) { + LOG(FATAL) << "No any nodes or elements"; + } + nodeDofs.initDofTable(nNodes()); elementDofs.initDofTable(nElements()); - for (uint32 el = 0; el < nElements(); el++) { + for (uint32 el = 0; el < nElements(); el++) { elements[el]->pre(); } @@ -354,7 +354,7 @@ void FEStorage::initDofs() { } // Total number of dofs (only registered by elements) - _nDofs = elementDofs.getNumberOfUsedDofs() + nodeDofs.getNumberOfUsedDofs(); + _nDofs = elementDofs.getNumberOfUsedDofs() + nodeDofs.getNumberOfUsedDofs(); CHECK(_nDofs); auto udofs = getUniqueNodeDofTypes(); @@ -381,11 +381,11 @@ void FEStorage::initDofs() { void FEStorage::assignEquationNumbers() { // _nUnknownDofs - number of Dof need to be found on every step - _nUnknownDofs = nDofs() - nConstrainedDofs(); + _nUnknownDofs = nDofs() - nConstrainedDofs(); CHECK(_nUnknownDofs); - uint32 next_eq_solve = nConstrainedDofs() + 1; - uint32 next_eq_const = 1; + uint32 next_eq_solve = nConstrainedDofs() + 1; + uint32 next_eq_const = 1; // TODO: we can try to do numbering in more efficient way if will loop over element nodes (as it // does in buildK() procedure) @@ -427,13 +427,13 @@ void FEStorage::assignEquationNumbers() { assert(next_eq_solve - 1 == nDofs() + nMpc()); - LOG(INFO) << "DoFs = " << nDofs() << ", constrained DoFs = " << nConstrainedDofs() << ", MPC eq. = " + LOG(INFO) << "DoFs = " << nDofs() << ", constrained DoFs = " << nConstrainedDofs() << ", MPC eq. = " << nMpc() << ", TOTAL eq. = " << nUnknownDofs() + nMpc(); } void FEStorage::initSolutionData () { TIMED_SCOPE(t, "initSolutionData"); - + // We need to know topology of the mesh in order to determine SparsityInfo for sparse matrices learnTopology(); @@ -449,14 +449,14 @@ void FEStorage::initSolutionData () { // | | | * | Us | = | Fs | + | Rs | // |KcsMPCc| KssMPCs | |-------| | | | | // | ^T | | | Ul | | Fl | | Rl | - // + // // 1. But really only // // | | | Us | | | // | KssMPCs | * |-------| = |RHS | // | | | Ul | | | // - // will be be solved by eq. solver. + // will be be solved by eq. solver. // // where: // @@ -480,7 +480,7 @@ void FEStorage::initSolutionData () { } - // initialize solutions vectors + // initialize solutions vectors vecU.reinit(nDofs()+nMpc()); if (transient) { vecDU.reinit(nDofs()+nMpc()); @@ -489,9 +489,9 @@ void FEStorage::initSolutionData () { vecR.reinit(nDofs()+nMpc()); vecF.reinit(nDofs()+nMpc()); - if (!material) { - LOG(WARNING) << "FEStorage::initializeSolutionData: material isn't defined"; - } + if (!material) { + LOG(WARNING) << "FEStorage::initializeSolutionData: material isn't defined"; + } // Need to restore non-zero entries in Sparse Matrices based on mesh topology and registered Dofs @@ -554,14 +554,14 @@ void FEStorage::printDofInfo(std::ostream& out) { for (uint32 en = 1; en <= nElements(); en++) { auto en_dofs = elementDofs.getEntityDofs(en); for (auto d1 = en_dofs.first; d1 != en_dofs.second; d1++) - out << "E" << en << ":" << Dof::dofType2label(d1->type) << " eq = " << d1->eqNumber + out << "E" << en << ":" << Dof::dofType2label(d1->type) << " eq = " << d1->eqNumber << " constrained = " << d1->isConstrained << endl; } for (uint32 nn = 1; nn <= nNodes(); nn++) { auto nn_dofs = nodeDofs.getEntityDofs(nn); for (auto d1 = nn_dofs.first; d1 != nn_dofs.second; d1++) - out << "N" << nn << ":" << Dof::dofType2label(d1->type) << " eq = " << d1->eqNumber + out << "N" << nn << ":" << Dof::dofType2label(d1->type) << " eq = " << d1->eqNumber << " constrained = " << d1->isConstrained << endl; } @@ -589,4 +589,4 @@ void FEStorage::learnTopology() { } -} // namespace nla3d +} // namespace nla3d diff --git a/src/lib/FEStorage.h b/src/lib/FEStorage.h index 79f23f1..fb386ea 100755 --- a/src/lib/FEStorage.h +++ b/src/lib/FEStorage.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once #include @@ -10,7 +10,7 @@ #include "math/BlockSparseMatrix.h" #include "FEComponent.h" #include "Mpc.h" - + namespace nla3d { class Element; @@ -19,7 +19,7 @@ class Dof; class ElementFactory; -// FEStorage - heart of the nla3d program, it contains all FE data: +// FEStorage - heart of the nla3d program, it contains all FE data: // * nodes array `nodes`, elements array `elements`, // * nodal DoFs array `elementDofs`, element DoFs array `nodeDofs`, // * solution matrices and vectors (`matK`, `vecU`, `vecF`, `vecR` and others), @@ -37,8 +37,8 @@ class ElementFactory; // should start from 1 and up to nElements()/nNodes(). class FEStorage { public: - FEStorage(); - ~FEStorage(); + FEStorage(); + ~FEStorage(); // A pointer to Material instance. nla3d supports only one material for FE model. External code // should create material instance and pass it to FEStorage. Then FEStorage takes a control on @@ -88,11 +88,11 @@ class FEStorage { void addValueR(uint32 eqi, double value); // fill with zeros the matrices of global system of equations - void zeroK(); - void zeroC(); - void zeroM(); + void zeroK(); + void zeroC(); + void zeroM(); // fill with zeros the vector vecF - void zeroF(); + void zeroF(); // get the K, C, M matrices (stiffness, damping, inertia) of the global system of equations. In current realization the // matrix is symmetric, but semi-positive defined, because of MPC equations (zeros on diagonals in @@ -115,25 +115,25 @@ class FEStorage { void assembleGlobalEqMatrices(); // getters to get numbers of different entities stored in FEStorage - uint32 nNodes(); - uint32 nElements(); - uint32 nDofs(); - uint32 nUnknownDofs(); - uint32 nConstrainedDofs(); - uint32 nMpc(); + uint32 nNodes(); + uint32 nElements(); + uint32 nDofs(); + uint32 nUnknownDofs(); + uint32 nConstrainedDofs(); + uint32 nMpc(); // if isTransient() == bool then FEStorage initialize matC, matM, vecDU, vecDDU along with matK, - // vecU. + // vecU. void setTransient(bool _transient); bool isTransient(); // Operations with DoFs // // Registation of DoFs is a key moment in nla3d. Every element (and other entities like MPC - // collections) should register its DoFs which the element will use. + // collections) should register its DoFs which the element will use. // For example: addNodeDof(32, Dof::UX) means that DoF UX for node 32 will be used in a global system of // equations. Unregistered DoFs will not participate in global equations system. This mechanism - // leads to ability to have different DoFs in different nodes. + // leads to ability to have different DoFs in different nodes. // NOTE: `node` > 0 void addNodeDof(uint32 node, std::initializer_list _dofs); // NOTE: `element` > 0 @@ -145,8 +145,8 @@ class FEStorage { bool isElementDofUsed(uint32 el, Dof::dofType dof); bool isNodeDofUsed(uint32 node, Dof::dofType dof); // get the number of equation in the global eq. system for particular DoF - uint32 getElementDofEqNumber(uint32 el, Dof::dofType dof); - uint32 getNodeDofEqNumber(uint32 node, Dof::dofType dof); + uint32 getElementDofEqNumber(uint32 el, Dof::dofType dof); + uint32 getNodeDofEqNumber(uint32 node, Dof::dofType dof); // get a pointer to Dof class for particular DoF Dof* getElementDof(uint32 el, Dof::dofType dof); Dof* getNodeDof(uint32 node, Dof::dofType dof); @@ -155,7 +155,7 @@ class FEStorage { double getElementDofSolution(uint32 el, Dof::dofType dof); // as far as constrained (fixed) DoFs have special treatment in FEStorage outer code (commonly // FESolver) should register all constrained DoFs by mean of these functions (before calling - // assignEquationNumbers()) + // assignEquationNumbers()) void setConstrainedNodeDof(uint32 node, Dof::dofType dtype); void setConstrainedElementDof(uint32 el, Dof::dofType dtype); // Get value of reaction for particular DoFs. Reactions are defined only for constrained DoFs. If @@ -163,27 +163,27 @@ class FEStorage { // returned. // NOTE: actually this functions just return vecR[eq-1] value, for eq <= nConstrainedDofs() it's // reaction load, but for eq > nConstrainedDofs() it's external DoF's concentrated load. - double getReaction(uint32 node, Dof::dofType dof); - double getReaction(uint32 eq); + double getReaction(uint32 node, Dof::dofType dof); + double getReaction(uint32 eq); // get an instance of Material class - Material* getMaterial(); + Material* getMaterial(); // get an instance of particular node // NOTE: `_nn` > 0 - Node& getNode(uint32 _nn); + Node& getNode(uint32 _nn); // function fills node_ptr with pointers to Node classes for element el. // Calling side should reserve a space for an array of pointers node_ptr. // size of node_ptr should be at least Element::getNNodes() // NOTE: `el` > 0 - void getElementNodes(uint32 el, Node** node_ptr); + void getElementNodes(uint32 el, Node** node_ptr); // The function return a spatial position of node `n` in either initial state (`deformed` = false) // or deformed state (`deformed` = true); // A calling side should take care about memory allocation for ptr. Size of ptr should be at least // 3 as the function always return 3-dimensional coordinates. - void getNodePosition(uint32 n, double* ptr, bool deformed = false); + void getNodePosition(uint32 n, double* ptr, bool deformed = false); // NOTE: `_en` > 0 - Element& getElement(uint32 _en); + Element& getElement(uint32 _en); // get a FEComponent instance by registration number // NOTE: `i` > -1 FEComponent* getFEComponent(size_t i); @@ -193,7 +193,7 @@ class FEStorage { // storing operations // // Add `mpc` to `mpcs` array. FEStorage will delete Mpc instances by itslef. - void addMpc(Mpc* mpc); + void addMpc(Mpc* mpc); // Add `mpcCol` to `mpcCollections` array. FEStorage will delete MpcCollection instances by itslef. void addMpcCollection(MpcCollection* mpcCol); // Add `comp` to `feComponents` array. FEStorage will delete FEComponent instances by itslef. @@ -203,10 +203,10 @@ class FEStorage { // instanses in `nodes` array void addNode(Node* node); // The function creates an vector `nodes` with dynamically allocated Node instances inside. - // Numbers of newly created elements pass back to the caller. - // NOTE: nodes will be created with default Node() constructor (node coordinates 0, 0, 0). - // NOTE: new nodes are concantenated to old nodes which already were in FEStorage - std::vector createNodes(uint32 nn); + // Numbers of newly created elements pass back to the caller. + // NOTE: nodes will be created with default Node() constructor (node coordinates 0, 0, 0). + // NOTE: new nodes are concantenated to old nodes which already were in FEStorage + std::vector createNodes(uint32 nn); // add an element to the element array `elements`. // TODO: this is inefficient functions because we can't ensure memory localization for all // elements instances in `elements` array @@ -214,13 +214,13 @@ class FEStorage { // The function creates a vector with dynamically allocated Element instances inside. Particular // realization of abstract Element class is chosen from elType variable by mean of ElementFactory // class (see elements/ElemenetFactory.h). Numbers of newly created elements pass back to the - // caller. + // caller. // NOTE: elements will be created with default constructor. User code should then fill node // numbers (by Element::getNodeNumber(..)) and other stuff related to particular realization of // Element class. std::vector createElements(uint32 en, ElementType elType); // template version of function described above. User code should provide example Element entity - // of particular class. + // of particular class. template std::vector createElements(uint32 _en, T example); @@ -228,7 +228,7 @@ class FEStorage { // // delete all FE model things: elements, nodes, MPCs, MPC Collections, FE components, topology // info. - void deleteMesh(); + void deleteMesh(); // delete nodes table: delete all dynamically allocated Node instances, // and clear vector of pointers `nodes`. void deleteNodes(); @@ -254,18 +254,18 @@ class FEStorage { // These functions perform all preparing steps before FESolver will solve the problem: // After this procedure one can't add more elements, nodes, boundary conditions, MPCs, .. // 1. Fill DoFs tables by calling Element->pre() and MpcCollection->pre() to register needed DoFs. - void initDofs(); + void initDofs(); // 2. Outer code should call FEStorage::setConstrained[Node/Element]Dof(..) to set the particular // DoFs to be constrained (fixed). // 3. Assign number for globals system equations by calling assignEquationNumbers(); void assignEquationNumbers(); // 4. Allocate all solution data structures: matK/C/M, vecU/DU/DDU/F/R - void initSolutionData(); + void initSolutionData(); // After global equations system is solved and vecU/DU/DDU/R is updated with appropriate values // FESolver should call this procedure to update element solution data // NOTE: actually Element::update() is called - void updateResults(); + void updateResults(); private: // fill `topology` data based on the current mesh (Element::nodes numbers) @@ -278,34 +278,34 @@ class FEStorage { void addEntryMPC(uint32 eq_num, uint32 eqj); // Total number of DoFs (registered by FEStorage::add[Node/Element]Dof(..)) - uint32 _nDofs = 0; + uint32 _nDofs = 0; // Number of constrained (fixed) DoFs values (registered by // FEStorage::setConstrained[Node/Element]Dof()) - uint32 _nConstrainedDofs = 0; + uint32 _nConstrainedDofs = 0; // _nUnknownDofs - number of Dof need to be found on every step (excluding Lagrangian lambdas) - uint32 _nUnknownDofs; + uint32 _nUnknownDofs; // Array of elements. Elements are created by the FEStorage::createElements(..) function. It's // important that elements in nla3d have consecutive numbering. That means that they have numbers // in [1; nElements()]. // NOTE: Elements are deleted by deleteElements() function. - std::vector elements; + std::vector elements; // Array of nodes. Nodes are created by FEStorage::createNodes(..). Nodes are created with default - // coordinates (0,0,0). Nodes have consecutive numbering. They have numbers in [1; nNodes()]. - std::vector nodes; + // coordinates (0,0,0). Nodes have consecutive numbering. They have numbers in [1; nNodes()]. + std::vector nodes; // List of MPC equations. List is populated by FEStorage::addMpc(..). An instance of Mpc class is // created outside of FEStorage class. But after addMpc(..) function FEStorage takes control on // the instance of Mpc. MPCs can be deleted by FEStorage::deleteMpcs(). Of course, they are - // deleted in ~FEStorage() destructor. - std::list mpcs; + // deleted in ~FEStorage() destructor. + std::list mpcs; // FE components - lists consist of numbers of entities. Here is the support of nodal components // and element components. This components can be used to apply BCs and/or MPCs. FE component is // attached to FE Storage by FEStorage::addFEComponent(..) function. All stored components can be // listed by FEStorage::listFEComponents(..). One can access to the particular component by - // methods FEStorage:: getFEComponent(..). + // methods FEStorage:: getFEComponent(..). std::vector feComponents; // An array of MpcCollection. Some number of MPCs can be grouped together because of they were @@ -322,7 +322,7 @@ class FEStorage { // only block(2) should be solved by eq. solver. // NOTE: Such separation is possible because FEStorage gives to constrained DoFs equation numbers // from 1 to nConstrainedDofs(). Unknown DoFs has numbers from nConstrainedDofs() + 1 to - // nConstrainedDofs() + nUnknownDofs(). + // nConstrainedDofs() + nUnknownDofs(). math::BlockSparseSymMatrix<2>* matK = nullptr; // matrices for transient analysis: @@ -337,7 +337,7 @@ class FEStorage { // this vector. // The whole vector could be divided into different parts: values for constrained DoFs in range // [0; nConstrainedDofs()), values for unknown DoFs [nConstrainedDofs(); nDofs()), values for - // Mpc's Lagrangian lambdas [nDofs(); nDofs() + nMpc()). + // Mpc's Lagrangian lambdas [nDofs(); nDofs() + nMpc()). math::dVec vecU; // Found first time derivatives for DoF's values. FESovler is responsible of filling this vector. @@ -350,7 +350,7 @@ class FEStorage { // this vector. // The whole vector could be divided into different parts: values of reactions for constrained DoFs in range // [0; nConstrainedDofs()), values of external loads for unknown DoFs [nConstrainedDofs(); nDofs()), for - // Mpc's Lagrangian lambdas [nDofs(); nDofs() + nMpc()) values should be alway remain zeros. + // Mpc's Lagrangian lambdas [nDofs(); nDofs() + nMpc()) values should be alway remain zeros. math::dVec vecR; // internal element loads. vecF has size nDofs + nMpc(). FEStorage::assembleGlobalEqMatrices() is @@ -359,7 +359,7 @@ class FEStorage { // * In nla3d every DoF is represented by `Dof` class. // * DoF can be nodal or for an element. - // * `elementDofs` and `nodeDofs` are DofCollection + // * `elementDofs` and `nodeDofs` are DofCollection // DofCollection keeps Dof objects and give it by (node, dof) pair // * Elements (or another entities) should register DoFs by using // FEStorage::add[Node/Element]Dof(node, dof) @@ -377,49 +377,49 @@ class FEStorage { inline void FEStorage::addValueK(uint32 nodei, Dof::dofType dofi, uint32 nodej, Dof::dofType dofj, double value) { - uint32 rowEq = getNodeDofEqNumber(nodei, dofi); - uint32 colEq = getNodeDofEqNumber(nodej, dofj); + uint32 rowEq = getNodeDofEqNumber(nodei, dofi); + uint32 colEq = getNodeDofEqNumber(nodej, dofj); addValueK(rowEq, colEq, value); } inline void FEStorage::addValueK(uint32 eqi, uint32 eqj, double value) { - // eqi - row equation + // eqi - row equation // eqj - column equation matK->addValue(eqi, eqj, value); } inline void FEStorage::addValueC(uint32 nodei, Dof::dofType dofi, uint32 nodej, Dof::dofType dofj, double value) { - uint32 rowEq = getNodeDofEqNumber(nodei, dofi); - uint32 colEq = getNodeDofEqNumber(nodej, dofj); + uint32 rowEq = getNodeDofEqNumber(nodei, dofi); + uint32 colEq = getNodeDofEqNumber(nodej, dofj); addValueC(rowEq, colEq, value); } inline void FEStorage::addValueC(uint32 eqi, uint32 eqj, double value) { - // eqi - row equation + // eqi - row equation // eqj - column equation matC->addValue(eqi, eqj, value); } inline void FEStorage::addValueM(uint32 nodei, Dof::dofType dofi, uint32 nodej, Dof::dofType dofj, double value) { - uint32 rowEq = getNodeDofEqNumber(nodei, dofi); - uint32 colEq = getNodeDofEqNumber(nodej, dofj); + uint32 rowEq = getNodeDofEqNumber(nodei, dofi); + uint32 colEq = getNodeDofEqNumber(nodej, dofj); addValueM(rowEq, colEq, value); } inline void FEStorage::addValueM(uint32 eqi, uint32 eqj, double value) { - // eqi - row equation + // eqi - row equation // eqj - column equation matM->addValue(eqi, eqj, value); } inline void FEStorage::addValueMPC(uint32 eq_num, uint32 nodej, Dof::dofType dofj, double coef) { - uint32 colEq = getNodeDofEqNumber(nodej, dofj); + uint32 colEq = getNodeDofEqNumber(nodej, dofj); addValueMPC(eq_num, colEq, coef); } @@ -430,47 +430,47 @@ inline void FEStorage::addValueMPC(uint32 eq_num, uint32 eqj, double coef) { inline void FEStorage::addValueF(uint32 nodei, Dof::dofType dofi, double value) { - uint32 rowEq = getNodeDofEqNumber(nodei, dofi); + uint32 rowEq = getNodeDofEqNumber(nodei, dofi); addValueF(rowEq, value); } inline void FEStorage::addValueF(uint32 eqi, double value) { assert(eqi > 0); - assert(eqi <= vecF.size()); - assert(eqi <= nDofs() + nMpc()); - vecF[eqi - 1] += value; + assert(eqi <= vecF.size()); + assert(eqi <= nDofs() + nMpc()); + vecF[eqi - 1] += value; } inline void FEStorage::addValueR(uint32 nodei, Dof::dofType dofi, double value) { - uint32 rowEq = getNodeDofEqNumber(nodei, dofi); + uint32 rowEq = getNodeDofEqNumber(nodei, dofi); addValueR(rowEq, value); } inline void FEStorage::addValueR(uint32 eqi, double value) { assert(eqi > 0); - assert(eqi <= vecR.size()); - assert(eqi <= nDofs() + nMpc()); - vecR[eqi - 1] += value; + assert(eqi <= vecR.size()); + assert(eqi <= nDofs() + nMpc()); + vecR[eqi - 1] += value; } inline void FEStorage::zeroK() { - assert(matK); + assert(matK); matK->zero(); } inline void FEStorage::zeroC() { - assert(matC); + assert(matC); matC->zero(); } inline void FEStorage::zeroM() { - assert(matM); + assert(matM); matM->zero(); } @@ -481,20 +481,20 @@ inline void FEStorage::zeroF() { inline math::BlockSparseSymMatrix<2>* FEStorage::getK() { - assert(matK->isCompressed()); - return matK; + assert(matK->isCompressed()); + return matK; } inline math::BlockSparseSymMatrix<2>* FEStorage::getC() { - assert(matC->isCompressed()); - return matC; + assert(matC->isCompressed()); + return matC; } inline math::BlockSparseSymMatrix<2>* FEStorage::getM() { - assert(matM->isCompressed()); - return matM; + assert(matM->isCompressed()); + return matM; } @@ -631,34 +631,34 @@ inline double FEStorage::getElementDofSolution (uint32 el, Dof::dofType dof) { inline Node& FEStorage::getNode(uint32 _nn) { - assert(_nn> 0 && _nn <= nNodes()); - return *nodes[_nn-1]; + assert(_nn> 0 && _nn <= nNodes()); + return *nodes[_nn-1]; } inline Element& FEStorage::getElement(uint32 _en) { - assert(_en <= nElements()); - return *(elements[_en-1]); + assert(_en <= nElements()); + return *(elements[_en-1]); } inline void FEStorage::addEntryK(uint32 eqi, uint32 eqj) { - // eqi - row equation + // eqi - row equation // eqj - column equation matK->addEntry(eqi, eqj); } inline void FEStorage::addEntryMPC(uint32 eq_num, uint32 eqj) { - assert(matK); - assert(eq_num > nConstrainedDofs() + nUnknownDofs()); + assert(matK); + assert(eq_num > nConstrainedDofs() + nUnknownDofs()); assert(eq_num <= nConstrainedDofs() + nUnknownDofs() + nMpc()); assert(eqj > 0 && eqj <= nConstrainedDofs() + nUnknownDofs()); matK->addEntry(eq_num, eqj); } -} // namespace nla3d +} // namespace nla3d // 'dirty' hack to avoid include loops (element-vs-festorage) #include "elements/element.h" @@ -689,4 +689,4 @@ std::vector FEStorage::createElements(uint32 _en, T example) { -} // namespace nla3d +} // namespace nla3d diff --git a/src/lib/Mpc.cpp b/src/lib/Mpc.cpp index 2d48436..f83522a 100644 --- a/src/lib/Mpc.cpp +++ b/src/lib/Mpc.cpp @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #include "Mpc.h" #include "FEStorage.h" @@ -76,7 +76,7 @@ void RigidBodyMpc::update () { // {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) + // 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); @@ -109,7 +109,7 @@ void RigidBodyMpc::update () { 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 + + C[i][j] = cos(thetaNorm) * solidmech::I[i][j] + c1 * theta0Mat_ij + c2 * theta0[i] * theta0[j]; } } diff --git a/src/lib/Mpc.h b/src/lib/Mpc.h index 22c4a5b..c5285a2 100644 --- a/src/lib/Mpc.h +++ b/src/lib/Mpc.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once #include diff --git a/src/lib/Node.cpp b/src/lib/Node.cpp index 3a8a673..30dc05c 100644 --- a/src/lib/Node.cpp +++ b/src/lib/Node.cpp @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #include "Node.h" diff --git a/src/lib/Node.h b/src/lib/Node.h index aba9372..e439193 100644 --- a/src/lib/Node.h +++ b/src/lib/Node.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once #include "sys.h" @@ -21,4 +21,4 @@ class Node { friend class Element; }; -} // namespace nla3d +} // namespace nla3d diff --git a/src/lib/PostProcessor.cpp b/src/lib/PostProcessor.cpp index 43c3477..ee5cf2c 100755 --- a/src/lib/PostProcessor.cpp +++ b/src/lib/PostProcessor.cpp @@ -1,22 +1,22 @@ // 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 +// https://github.com/dmitryikh/nla3d #include "PostProcessor.h" namespace nla3d { PostProcessor::PostProcessor(FEStorage *st) { - assert(st); - storage = st; - active = false; - failed = false; + 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); + 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 +} // namespace nla3d diff --git a/src/lib/PostProcessor.h b/src/lib/PostProcessor.h index 3ad5b30..b8a46c8 100755 --- a/src/lib/PostProcessor.h +++ b/src/lib/PostProcessor.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once #include "sys.h" @@ -14,39 +14,39 @@ 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," + 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; - } + return; + } + active = act; + } + bool getActive () { + return active; + } - uint16 getProc_num () { - return nPost_proc; - } + uint16 getProc_num () { + return nPost_proc; + } - friend class FEStorage; - friend class FESolver; + friend class FEStorage; + friend class FESolver; protected: - FEStorage *storage; - uint16 nPost_proc; - std::string name; - bool active; - bool failed; + FEStorage *storage; + uint16 nPost_proc; + std::string name; + bool active; + bool failed; }; } // namespace nla3d diff --git a/src/lib/ReactionProcessor.cpp b/src/lib/ReactionProcessor.cpp index 72ab1b7..0d974e9 100755 --- a/src/lib/ReactionProcessor.cpp +++ b/src/lib/ReactionProcessor.cpp @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #include "ReactionProcessor.h" #include "Node.h" @@ -8,16 +8,16 @@ namespace nla3d { ReactionProcessor::ReactionProcessor(FEStorage *st) : PostProcessor(st) { - name ="ReactionProcessor"; + name ="ReactionProcessor"; } ReactionProcessor::ReactionProcessor(FEStorage *st, std::string _filename) : PostProcessor(st) { - name ="ReactionProcessor"; + name ="ReactionProcessor"; filename = _filename; } void ReactionProcessor::pre() { - if (nodes.size() == 0) { + if (nodes.size() == 0) { LOG(WARNING) << "Can't work. No nodes. Processor name = " << name; return; } @@ -31,11 +31,11 @@ void ReactionProcessor::pre() { Dof::dofType t = static_cast (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; @@ -92,7 +92,7 @@ void ReactionProcessor::process (uint16 curLoadstep) { std::vector reactions; reactions.assign(dofs.size(), 0.0); - for (uint32 n = 0; n < nodes.size(); n++) { + for (uint32 n = 0; n < nodes.size(); n++) { for (uint16 d = 0; d < dofs.size(); d++) { reactions[d] += storage->getReaction(nodes[n],dofs[d]); } diff --git a/src/lib/ReactionProcessor.h b/src/lib/ReactionProcessor.h index 94b1a3c..2415e1a 100755 --- a/src/lib/ReactionProcessor.h +++ b/src/lib/ReactionProcessor.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once #include "sys.h" @@ -11,21 +11,21 @@ namespace nla3d { class ReactionProcessor : public PostProcessor { public: - ReactionProcessor(FEStorage *st); - ReactionProcessor(FEStorage *st, std::string _filename); - virtual ~ReactionProcessor() { }; + ReactionProcessor(FEStorage *st); + ReactionProcessor(FEStorage *st, std::string _filename); + virtual ~ReactionProcessor() { }; - virtual void pre (); - virtual void process (uint16 curLoadstep); - virtual void post (uint16 curLoadstep); + virtual void pre (); + virtual void process (uint16 curLoadstep); + virtual void post (uint16 curLoadstep); std::vector getReactions (Dof::dofType dof); - - std::vector nodes; - std::vector dofs; + + std::vector nodes; + std::vector dofs; protected: std::string filename; - std::vector > sumOfDofsReactions; + std::vector > sumOfDofsReactions; }; } // namespace nla3d diff --git a/src/lib/VtkProcessor.cpp b/src/lib/VtkProcessor.cpp index 4dacf72..9fd47a9 100755 --- a/src/lib/VtkProcessor.cpp +++ b/src/lib/VtkProcessor.cpp @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #include "VtkProcessor.h" #include "elements/element.h" @@ -9,8 +9,8 @@ namespace nla3d { using namespace math; VtkProcessor::VtkProcessor(FEStorage *st, std::string _fileName) : PostProcessor(st) { - name = "VtkProcessor"; - file_name = _fileName; + name = "VtkProcessor"; + file_name = _fileName; } VtkProcessor::~VtkProcessor() { @@ -46,23 +46,23 @@ void VtkProcessor::pre(){ } // write zero VTK file (before solution) - std::string cur_fn = file_name + "0" + ".vtk"; - std::ofstream file(cur_fn.c_str(), std::ios::trunc); - write_header(file); - write_geometry(file); - write_point_data(file); - write_cell_data(file); - file.close(); + std::string cur_fn = file_name + "0" + ".vtk"; + std::ofstream file(cur_fn.c_str(), std::ios::trunc); + write_header(file); + write_geometry(file); + write_point_data(file); + write_cell_data(file); + file.close(); } void VtkProcessor::process (uint16 curLoadstep) { - std::string cur_fn = file_name + toStr(curLoadstep) + ".vtk"; - std::ofstream file(cur_fn.c_str(), std::ios::trunc); - write_header(file); - write_geometry(file); - write_point_data(file); - write_cell_data(file); - file.close(); + std::string cur_fn = file_name + toStr(curLoadstep) + ".vtk"; + std::ofstream file(cur_fn.c_str(), std::ios::trunc); + write_header(file); + write_geometry(file); + write_point_data(file); + write_cell_data(file); + file.close(); } void VtkProcessor::post (uint16 curLoadstep) { @@ -108,44 +108,44 @@ bool VtkProcessor::registerResults(std::initializer_list queries) { void VtkProcessor::write_header (std::ofstream &file) { - file << "# vtk DataFile Version 2.0" << std::endl; - file << "This file is generated by nla3D program" << std::endl << "ASCII" << std::endl; - file << "DATASET UNSTRUCTURED_GRID" << std::endl; + file << "# vtk DataFile Version 2.0" << std::endl; + file << "This file is generated by nla3D program" << std::endl << "ASCII" << std::endl; + file << "DATASET UNSTRUCTURED_GRID" << std::endl; } void VtkProcessor::write_geometry(std::ofstream &file, bool def) { - Vec<3> xi; - file << "POINTS " << storage->nNodes() << " float" << std::endl; - for (uint32 i=1; i <= storage->nNodes(); i++) - { - storage->getNodePosition(i, xi.ptr(), def); - file << xi << std::endl; - } - /* - CELLS en en*9 - 4 i j k l - - CELL_TYPES en - 11 - 11 + Vec<3> xi; + file << "POINTS " << storage->nNodes() << " float" << std::endl; + for (uint32 i=1; i <= storage->nNodes(); i++) + { + storage->getNodePosition(i, xi.ptr(), def); + file << xi << std::endl; + } + /* + CELLS en en*9 + 4 i j k l + + CELL_TYPES en + 11 + 11 */ // calculate overall size of element data uint32 data_size = 0; - for (uint32 i=1; i <= storage->nElements(); i++) + for (uint32 i=1; i <= storage->nElements(); i++) data_size += storage->getElement(i).getNNodes() + 1; - file << "CELLS " << storage->nElements() << " " << data_size << std::endl; - for (uint32 i=1; i <= storage->nElements(); i++) - { + file << "CELLS " << storage->nElements() << " " << data_size << std::endl; + for (uint32 i=1; i <= storage->nElements(); i++) + { uint16 nodesNum = storage->getElement(i).getNNodes(); - file << nodesNum; - for (uint16 j=0; j < nodesNum; j++) - file << " " << storage->getElement(i).getNodeNumber(j)-1; - file << std::endl; - } - file << "CELL_TYPES " << storage->nElements() << std::endl; - for (uint32 i=1; i <= storage->nElements(); i++) { + file << nodesNum; + for (uint16 j=0; j < nodesNum; j++) + file << " " << storage->getElement(i).getNodeNumber(j)-1; + file << std::endl; + } + file << "CELL_TYPES " << storage->nElements() << std::endl; + for (uint32 i=1; i <= storage->nElements(); i++) { ElementShape eltype = storage->getElement(i).getShape(); if (eltype == ElementShape::QUAD) { file << "9" << std::endl; //VTK_QUAD, see VTK file formats @@ -207,8 +207,8 @@ void VtkProcessor::write_point_data(std::ofstream &file) { } //Write cell data averaging from all integration points. -//Use global coordinate system. all futher transformations -//should be done on paraview side. +//Use global coordinate system. all futher transformations +//should be done on paraview side. void VtkProcessor::write_cell_data(std::ofstream &file) { size_t en = storage->nElements(); std::vector > dataTensor; @@ -221,7 +221,7 @@ void VtkProcessor::write_cell_data(std::ofstream &file) { return; } - file << "CELL_DATA " << en << std::endl; + file << "CELL_DATA " << en << std::endl; // write element DoFs for (auto v : elementDofTypes) { diff --git a/src/lib/VtkProcessor.h b/src/lib/VtkProcessor.h index 84ce7e7..939db4d 100755 --- a/src/lib/VtkProcessor.h +++ b/src/lib/VtkProcessor.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once #include @@ -25,11 +25,11 @@ namespace nla3d { 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); + 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 @@ -43,14 +43,14 @@ class VtkProcessor : public PostProcessor { bool registerResults(std::initializer_list queries); private: - std::string file_name; + std::string file_name; - void write_header(std::ofstream &file); + void write_header(std::ofstream &file); // 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_point_data(std::ofstream &file); - void write_cell_data(std::ofstream &file); + void write_geometry(std::ofstream &file, bool def = false); + void write_point_data(std::ofstream &file); + void write_cell_data(std::ofstream &file); void writeScalar(std::ofstream &file, const char* name, std::vector& data); void writeTensor(std::ofstream &file, const char* name, std::vector >& data); diff --git a/src/lib/elements/ElementFactory.cpp b/src/lib/elements/ElementFactory.cpp index c6afb8f..1cbde90 100644 --- a/src/lib/elements/ElementFactory.cpp +++ b/src/lib/elements/ElementFactory.cpp @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #include "elements/ElementFactory.h" #include "elements/PLANE41.h" diff --git a/src/lib/elements/ElementFactory.h b/src/lib/elements/ElementFactory.h index 31e5c12..d13fe64 100644 --- a/src/lib/elements/ElementFactory.h +++ b/src/lib/elements/ElementFactory.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once #include "sys.h" @@ -9,10 +9,9 @@ namespace nla3d { class Element; class ElementFactory { - public: - - static ElementType elName2elType (std::string elName); - static void createElements (ElementType elId, const uint32 n, std::vector& ptr); +public: + static ElementType elName2elType (std::string elName); + static void createElements (ElementType elId, const uint32 n, std::vector& ptr); }; } // namespace nla3d diff --git a/src/lib/elements/PLANE41.cpp b/src/lib/elements/PLANE41.cpp index fe31d30..b2bb5e8 100755 --- a/src/lib/elements/PLANE41.cpp +++ b/src/lib/elements/PLANE41.cpp @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #include "elements/PLANE41.h" @@ -65,11 +65,11 @@ void ElementPLANE41::buildK() { Mat<3,8> matB = make_B(np); //матрица S для матричного умножения - Mat<4,4> matS = Mat<4,4>(S[np][0],S[np][2], 0.0, 0.0, + Mat<4,4> matS = Mat<4,4>(S[np][0],S[np][2], 0.0, 0.0, S[np][2],S[np][1], 0.0, 0.0, 0.0, 0.0, S[np][0],S[np][2], 0.0, 0.0, S[np][2],S[np][1]); - //матрица Омега.используется для составления + //матрица Омега.используется для составления //матр. накопленных линейных деформаций к текущему шагу Mat<3,4> matO = Mat<3,4>(O[np][0], 0.0, O[np][2], 0.0, 0.0, O[np][1], 0.0, O[np][3], @@ -85,7 +85,7 @@ void ElementPLANE41::buildK() { Kpp -= 1.0/k*dWt; }// loop over intergration points - + //сборка в одну матрицу for (uint16 i=0; i < 8; i++) @@ -152,7 +152,7 @@ void ElementPLANE41::update() { bool ElementPLANE41::getScalar(double* scalar, scalarQuery query, uint16 gp, const double scale) { //see queries in query.h - //gp - needed gauss point + //gp - needed gauss point assert(scalar != nullptr); if (gp == GP_MEAN) { //need to average result over the element @@ -254,8 +254,8 @@ bool ElementPLANE41::getTensor(math::MatSym<3>* tensor, tensorQuery query, uint J = matF.data[0][0]*(matF.data[1][1]*matF.data[2][2]-matF.data[1][2]*matF.data[2][1])- matF.data[0][1]*(matF.data[1][0]*matF.data[2][2]-matF.data[1][2]*matF.data[2][0])+ matF.data[0][2]*(matF.data[1][0]*matF.data[2][1]-matF.data[1][1]*matF.data[2][0]); - //In order to complete matS (3x3 symmetric matrix, PK2 tensor) we need - //to know S33 component: + //In order to complete matS (3x3 symmetric matrix, PK2 tensor) we need + //to know S33 component: //1) One solution is to calculate S33 on every solution step //and store it in S[np] vector. //2) Second solution is to resotre S33 right here. diff --git a/src/lib/elements/PLANE41.h b/src/lib/elements/PLANE41.h index 9639b52..840e84c 100755 --- a/src/lib/elements/PLANE41.h +++ b/src/lib/elements/PLANE41.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once #include "elements/element.h" @@ -96,4 +96,4 @@ void ElementPLANE41::assembleK(const math::Mat &Ke, con } } -} // namespace nla3d +} // namespace nla3d diff --git a/src/lib/elements/QUADTH.cpp b/src/lib/elements/QUADTH.cpp index 56da99b..89f994b 100644 --- a/src/lib/elements/QUADTH.cpp +++ b/src/lib/elements/QUADTH.cpp @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #include "elements/QUADTH.h" @@ -59,7 +59,7 @@ void ElementQUADTH::buildC() { dWt = intWeight(np); Vec<4> ff = formFunc(np); for (uint16 i = 0; i < 4; i++) - for (uint16 j = i; j < 4; j++) + for (uint16 j = i; j < 4; j++) Ce.comp(i, j) += ff[i] * ff[j] * rho * c * dWt; }// loop over integration points diff --git a/src/lib/elements/QUADTH.h b/src/lib/elements/QUADTH.h index ce26d5a..0c1e487 100644 --- a/src/lib/elements/QUADTH.h +++ b/src/lib/elements/QUADTH.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once #include "elements/element.h" diff --git a/src/lib/elements/SOLID81.cpp b/src/lib/elements/SOLID81.cpp index a3130f4..1b2ee2d 100644 --- a/src/lib/elements/SOLID81.cpp +++ b/src/lib/elements/SOLID81.cpp @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #include "elements/SOLID81.h" @@ -192,7 +192,7 @@ void ElementSOLID81::make_Omega (uint16 np, Mat<6,9> &B) { bool ElementSOLID81::getScalar(double* scalar, scalarQuery query, uint16 gp, const double scale) { //see queries in query.h - //gp - needed gauss point + //gp - needed gauss point if (gp == GP_MEAN) { //need to average result over the element double dWtSum = volume(); double dWt; @@ -288,7 +288,7 @@ bool ElementSOLID81::getTensor(math::MatSym<3>* tensor, tensorQuery query, uint switch (query) { case tensorQuery::COUCHY: - //matF^T + //matF^T matF.data[0][0] = 1+O[gp][0]; matF.data[1][0] = O[gp][1]; matF.data[2][0] = O[gp][2]; @@ -303,7 +303,7 @@ bool ElementSOLID81::getTensor(math::MatSym<3>* tensor, tensorQuery query, uint //deviatoric part of S: Sd = S[gp] //hydrostatic part of S: Sp = p * J * C^(-1) - solidmech::invC_C (C[gp].ptr(), J, cInv); + solidmech::invC_C (C[gp].ptr(), J, cInv); pe = storage->getElementDofSolution(getElNum(), Dof::HYDRO_PRESSURE); // TODO: it seems that S[gp] contains deviatoric + pressure already.. // we dont need to sum it again @@ -316,7 +316,7 @@ bool ElementSOLID81::getTensor(math::MatSym<3>* tensor, tensorQuery query, uint case tensorQuery::PK2: // hydrostatic part of S: Sp = p * J * C^(-1) J = solidmech::J_C(C[gp].ptr()); - solidmech::invC_C (C[gp].ptr(), J, cInv); + solidmech::invC_C (C[gp].ptr(), J, cInv); pe = storage->getElementDofSolution(getElNum(), Dof::HYDRO_PRESSURE); for (uint16 i = 0; i < 6; i++) { tensor->data[i] += (S[gp][i] + pe * J * cInv[i]) * scale; diff --git a/src/lib/elements/SOLID81.h b/src/lib/elements/SOLID81.h index 5407773..921f1eb 100644 --- a/src/lib/elements/SOLID81.h +++ b/src/lib/elements/SOLID81.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once #include "elements/element.h" @@ -26,7 +26,7 @@ class ElementSOLID81 : public ElementIsoParamHEXAHEDRON { void buildK(); void update(); - void make_B_L (uint16 nPoint, math::Mat<6,24> &B); //функция создает линейную матрицу [B] + void make_B_L (uint16 nPoint, math::Mat<6,24> &B); //функция создает линейную матрицу [B] void make_B_NL (uint16 nPoint, math::Mat<9,24> &B); //функция создает линейную матрицу [Bomega] void make_S (uint16 nPoint, math::MatSym<9> &B); void make_Omega (uint16 nPoint, math::Mat<6,9> &B); @@ -41,7 +41,7 @@ class ElementSOLID81 : public ElementIsoParamHEXAHEDRON { std::vector > S; //S[номер т. интегр.][номер напряжения] - напряжения Пиолы-Кирхгоффа //C[M_XX], C[M_XY], C[M_XZ], C[M_YY], C[M_YZ], C[M_ZZ] std::vector > C; //C[номер т. интегр.][номер деформ.] - компоненты тензора меры деформации - // O[0]-dU/dx O[1]-dU/dy O[2]-dU/dz O[3]-dV/dx O[4]-dV/dy O[5]-dV/dz O[6]-dW/dx O[7]-dW/dy O[8]-dW/dz + // O[0]-dU/dx O[1]-dU/dy O[2]-dU/dz O[3]-dV/dx O[4]-dV/dy O[5]-dV/dz O[6]-dW/dx O[7]-dW/dy O[8]-dW/dz std::vector > O; //S[номер т. интегр.][номер омеги] template @@ -54,41 +54,41 @@ class ElementSOLID81 : public ElementIsoParamHEXAHEDRON { template void ElementSOLID81::assemble3(math::MatSym &Kuu, math::Vec &Kup, double Kpp, math::Vec &Fu, double Fp) { const uint16 dim = 3; - assert (getNNodes() * dim == dimM); - double *Kuu_p = Kuu.ptr(); - double *Kup_p = Kup.ptr(); - double *Fu_p = Fu.ptr(); + assert (getNNodes() * dim == dimM); + double *Kuu_p = Kuu.ptr(); + double *Kup_p = Kup.ptr(); + double *Fu_p = Fu.ptr(); Dof::dofType dofVec[] = {Dof::UX, Dof::UY, Dof::UZ}; - for (uint16 i=0; i < getNNodes(); i++) - for (uint16 di=0; di < dim; di++) - for (uint16 j=i; j < getNNodes(); j++) - - for (uint16 dj=0; dj < dim; dj++) - { - if ((i==j) && (djaddValueK(nodes[i],dofVec[di],nodes[j], dofVec[dj], *Kuu_p); - Kuu_p++; - } - } - //upper diagonal process for nodes-el dofs + for (uint16 i=0; i < getNNodes(); i++) + for (uint16 di=0; di < dim; di++) + for (uint16 j=i; j < getNNodes(); j++) + + for (uint16 dj=0; dj < dim; dj++) + { + if ((i==j) && (djaddValueK(nodes[i],dofVec[di],nodes[j], dofVec[dj], *Kuu_p); + Kuu_p++; + } + } + //upper diagonal process for nodes-el dofs uint32 elEq = storage->getElementDofEqNumber(getElNum(), Dof::HYDRO_PRESSURE); - for (uint16 i=0; i < getNNodes(); i++) { - for(uint16 di=0; di < dim; di++) { + for (uint16 i=0; i < getNNodes(); i++) { + for(uint16 di=0; di < dim; di++) { uint32 rowEq = storage->getNodeDofEqNumber(nodes[i], dofVec[di]); storage->addValueK(rowEq, elEq, *Kup_p); Kup_p++; } } - //upper diagonal process for el-el dofs + //upper diagonal process for el-el dofs storage->addValueK(elEq, elEq, Kpp); - for (uint16 i=0; i < getNNodes(); i++) { - for (uint16 di=0; di < dim; di++) { - storage->addValueF(nodes[i], dofVec[di], *Fu_p); - Fu_p++; - } + for (uint16 i=0; i < getNNodes(); i++) { + for (uint16 di=0; di < dim; di++) { + storage->addValueF(nodes[i], dofVec[di], *Fu_p); + Fu_p++; + } } storage->addValueF(elEq, Fp); } diff --git a/src/lib/elements/TETRA0.cpp b/src/lib/elements/TETRA0.cpp index 924fc02..a8053fc 100644 --- a/src/lib/elements/TETRA0.cpp +++ b/src/lib/elements/TETRA0.cpp @@ -1,3 +1,7 @@ +// 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/TETRA0.h" using namespace std; @@ -39,10 +43,10 @@ void ElementTETRA0::buildK() { // fill here matC makeC(matC); // fill here matB - makeB(matB); + makeB(matB); math::matBTDBprod(matB, matC, vol, matKe); - // start assemble procedure. Here we should provide element stiffness matrix and an order of + // start assemble procedure. Here we should provide element stiffness matrix and an order of // nodal DoFs in the matrix. assembleK(matKe, {Dof::UX, Dof::UY, Dof::UZ}); } @@ -68,7 +72,7 @@ void ElementTETRA0::update () { U[i*3 + 1] = storage->getNodeDofSolution(getNodeNumber(i), Dof::UY); U[i*3 + 2] = storage->getNodeDofSolution(getNodeNumber(i), Dof::UZ); } - + // restore strains math::matBVprod(matB, U, -1.0, strains); math::matBVprod(matC, strains, 1.0, stress); @@ -168,7 +172,7 @@ void ElementTETRA0::makeC (math::MatSym<6> &C) { int ElementTETRA0::permute(int i){ if (i > 3) return i -4; - else + else return i; } @@ -200,7 +204,7 @@ bool ElementTETRA0::getTensor(math::MatSym<3>* tensor, tensorQuery query, uint1 return true; } - + return false; } } //namespace nla3d diff --git a/src/lib/elements/TETRA0.h b/src/lib/elements/TETRA0.h index 339a82d..e1f23e8 100644 --- a/src/lib/elements/TETRA0.h +++ b/src/lib/elements/TETRA0.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once #include "elements/element.h" diff --git a/src/lib/elements/TRIANGLE4.cpp b/src/lib/elements/TRIANGLE4.cpp index 41e2d7a..b795913 100644 --- a/src/lib/elements/TRIANGLE4.cpp +++ b/src/lib/elements/TRIANGLE4.cpp @@ -1,5 +1,8 @@ -#include "elements/TRIANGLE4.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 +#include "elements/TRIANGLE4.h" #include using namespace std; @@ -30,7 +33,7 @@ void ElementTRIANGLE4::buildK() { // matC is 3d elastic matrix math::MatSym<3> matC; matC.zero(); - //only for area + //only for area Eigen::MatrixXd matS(3,3); matS.setZero(); matS << 1. , storage->getNode(getNodeNumber(0)).pos[0] , storage->getNode(getNodeNumber(0)).pos[1] , @@ -41,9 +44,9 @@ void ElementTRIANGLE4::buildK() { makeC(matC); // fill here matB makeB(matB); - + math::matBTDBprod(matB, matC, area, matKe); - // start assemble procedure. Here we should provide element stiffness matrix and an order of + // start assemble procedure. Here we should provide element stiffness matrix and an order of // nodal DoFs in the matrix. assembleK(matKe, {Dof::UX, Dof::UY}); } @@ -113,9 +116,9 @@ void ElementTRIANGLE4::makeC (math::MatSym<3> &C) { if (state == PlaneState::Stress){ const double A = E*(1.-my*my); C.comp(0,0) = 1.*A; - C.comp(0,1) = my*A; + C.comp(0,1) = my*A; C.comp(1,1) = 1.*A; - C.comp(2,2) = (1.-my)/2.*A; + C.comp(2,2) = (1.-my)/2.*A; } } diff --git a/src/lib/elements/TRIANGLE4.h b/src/lib/elements/TRIANGLE4.h index e275941..58f4a01 100644 --- a/src/lib/elements/TRIANGLE4.h +++ b/src/lib/elements/TRIANGLE4.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once #include "elements/element.h" diff --git a/src/lib/elements/TRUSS3.cpp b/src/lib/elements/TRUSS3.cpp index 9bde419..99350f0 100644 --- a/src/lib/elements/TRUSS3.cpp +++ b/src/lib/elements/TRUSS3.cpp @@ -1,3 +1,7 @@ +// 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/TRUSS3.h" namespace nla3d { @@ -28,7 +32,7 @@ void ElementTRUSS3::buildK() { // here we use pointer to FEStorage class where all FE information is stored (from node and // element tables to solution results). getNode(..).pos used to obtain an class instance for - // particular node an get its position which is Vec<3> data type. + // particular node an get its position which is Vec<3> data type. math::Vec<3> deltaPos = storage->getNode(getNodeNumber(1)).pos - storage->getNode(getNodeNumber(0)).pos; // as long as we need 1.0/length instead of length itself it's useful to store this value in the // variable inv_length. Such kind of code optimization is quite useful in massively computational @@ -45,7 +49,7 @@ void ElementTRUSS3::buildK() { T(1,4) = T(0,1); T(1,5) = T(0,2); - // This is just another way in Eigen to initialize {{1.0, -1.0}, {-1.0, 1.0}} matrix. + // This is just another way in Eigen to initialize {{1.0, -1.0}, {-1.0, 1.0}} matrix. K << 1.0, -1.0, -1.0, 1.0; K *= A*E*inv_length; @@ -54,7 +58,7 @@ void ElementTRUSS3::buildK() { // Ke.triangularView is used. Ke.triangularView() = T.transpose() * K * T; - // start assemble procedure. Here we should provide element stiffness matrix and an order of + // start assemble procedure. Here we should provide element stiffness matrix and an order of // nodal DoFs in the matrix. assembleK(Ke, {Dof::UX, Dof::UY, Dof::UZ}); } diff --git a/src/lib/elements/TRUSS3.h b/src/lib/elements/TRUSS3.h index 2557554..4ee9cf4 100644 --- a/src/lib/elements/TRUSS3.h +++ b/src/lib/elements/TRUSS3.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once #include "elements/element.h" @@ -10,7 +10,7 @@ namespace nla3d { // ElementTRUSS3 (TRUSS3) is a simplest realisation of 3D truss finite element (FE). Its primiral goal is to // show how to implement new FE into nla3d library. One can find a lot of information about how to // build 3D truss element stiffnes matrix in the internet. The author was looking at the material on -// this website: +// this website: // http://what-when-how.com/the-finite-element-method/fem-for-trusses-finite-element-method-part-1/ // New FE formulation is encapsulated into a class derived from base Element class. FEs have a lot diff --git a/src/lib/elements/element.cpp b/src/lib/elements/element.cpp index b195a4d..9ee211a 100755 --- a/src/lib/elements/element.cpp +++ b/src/lib/elements/element.cpp @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #include "elements/element.h" @@ -50,7 +50,7 @@ void Element::assembleK(Eigen::Ref Ke, if ((i==j) && (djaddValueK(nodes[i], nodeDof[di], nodes[j], nodeDof[dj], + storage->addValueK(nodes[i], nodeDof[di], nodes[j], nodeDof[dj], Ke.selfadjointView()(i*dim+di, j*dim +dj)); } } diff --git a/src/lib/elements/element.h b/src/lib/elements/element.h index 2f3c2b2..45ad8e6 100755 --- a/src/lib/elements/element.h +++ b/src/lib/elements/element.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once #include "sys.h" @@ -81,7 +81,7 @@ static const uint16 _shape_nnodes[] = { // (methods getScalar(...), getVector(...), getTensor(...)), for update element state after solution // iteration (method update()). class Element { - public: + public: Element (); virtual ~Element(); @@ -89,7 +89,7 @@ class Element { uint32 getElNum(); // return number of nodes for the element uint16 getNNodes(); - // return number of dimensions (0D, 1D, 2D, 3D) occupied by element shape. + // return number of dimensions (0D, 1D, 2D, 3D) occupied by element shape. uint16 getDim(); // return element shape ElementShape getShape(); @@ -378,7 +378,7 @@ inline uint16 Element::getIntegrationOrder() { inline void Element::setIntegrationOrder(uint16 _nint) { - intOrder = _nint; + intOrder = _nint; } diff --git a/src/lib/elements/isoparametric.cpp b/src/lib/elements/isoparametric.cpp index 5c527e3..6f682cb 100644 --- a/src/lib/elements/isoparametric.cpp +++ b/src/lib/elements/isoparametric.cpp @@ -1,7 +1,6 @@ - // 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 +// https://github.com/dmitryikh/nla3d #include "elements/element.h" #include "isoparametric.h" @@ -91,7 +90,7 @@ void ElementIsoParamQUAD::makeJacob() { } det[np] = J.det(); // determinant of Jacob matrix - // check for geometry form error + // check for geometry form error LOG_IF(det[np] < 1.0e-20, ERROR) << "Determinant is too small (" << det[np] << ") in element = " << elNum; // обращение матрицы Якоби inv_det = 1.0/det[np]; @@ -173,7 +172,7 @@ void ElementIsoParamHEXAHEDRON::makeJacob() { } det[np] = J.det(); // determinant of Jacob matrix - // check for geometry form error + // check for geometry form error LOG_IF(det[np] < 1.0e-20, ERROR) << "Determinant is too small (" << det[np] << ") in element = " << elNum; // обращение матрицы Якоби inv_det = 1.0/det[np]; diff --git a/src/lib/elements/isoparametric.h b/src/lib/elements/isoparametric.h index 4de2e48..f2137da 100644 --- a/src/lib/elements/isoparametric.h +++ b/src/lib/elements/isoparametric.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once #include "sys.h" @@ -221,7 +221,7 @@ class ElementIsoParamQUAD : public ElementQUAD { double sideDet[4]; //Jacobian for side integration // function to calculate all staff for isoparametric FE - void makeJacob(); + void makeJacob(); double intWeight(uint16 np); double volume(); void np2rst(uint16 np, double *xi); //by number of gauss point find local coordinates @@ -244,7 +244,7 @@ class ElementIsoParamHEXAHEDRON : public ElementHEXAHEDRON { std::vector det; //Jacobian // function to calculate all staff for isoparametric FE - void makeJacob(); + void makeJacob(); double intWeight(uint16 np); double volume(); void np2rst(uint16 np, double *xi); //by number of gauss point find local coordinates diff --git a/src/lib/elements/query.h b/src/lib/elements/query.h index 038102c..4f9852d 100644 --- a/src/lib/elements/query.h +++ b/src/lib/elements/query.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once @@ -35,15 +35,15 @@ static_assert((int)vectorQuery::LAST == sizeof(vectorQueryLabels)/sizeof(vectorQ enum class tensorQuery { - UNDEF, + UNDEF, // usual stress tensor - COUCHY, + COUCHY, // second Piola-Kirchgoff stress tensor (symmetric 3x3) - PK2, + PK2, // Lagrange deformations - E, - // C = F^T F - C, + E, + // C = F^T F + C, LAST }; diff --git a/src/lib/materials/MaterialFactory.cpp b/src/lib/materials/MaterialFactory.cpp index 0e2284c..1cad5c2 100644 --- a/src/lib/materials/MaterialFactory.cpp +++ b/src/lib/materials/MaterialFactory.cpp @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #include "materials/MaterialFactory.h" diff --git a/src/lib/materials/MaterialFactory.h b/src/lib/materials/MaterialFactory.h index f0d1bef..c0088b6 100644 --- a/src/lib/materials/MaterialFactory.h +++ b/src/lib/materials/MaterialFactory.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once #include "materials/materials_hyperelastic.h" @@ -19,8 +19,8 @@ class MaterialFactory { static const char* const matModelLabels[]; - static matId matName2matId (std::string matName); - static Material* createMaterial (std::string matName); + static matId matName2matId (std::string matName); + static Material* createMaterial (std::string matName); }; } // namespace nla3d diff --git a/src/lib/materials/material.cpp b/src/lib/materials/material.cpp index bf31414..5333ce9 100755 --- a/src/lib/materials/material.cpp +++ b/src/lib/materials/material.cpp @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #include "materials/material.h" @@ -14,17 +14,17 @@ const double Material::I[6] = {1.0,0.0,0.0,1.0,0.0,1.0}; Material::Material (uint16 num_c) { code = 0; - numC = num_c; - MC = new double[num_c]; + numC = num_c; + MC = new double[num_c]; } std::string Material::toString() { std::string str; - str += getName(); - for (size_t i=0; i < getNumC(); i++) { - str += " " + MC_names[i] + " = " + toStr(MC[i]); + str += getName(); + for (size_t i=0; i < getNumC(); i++) { + str += " " + MC_names[i] + " = " + toStr(MC[i]); } - return str; + return str; } double& Material::Ci (const std::string& nameConst) { @@ -33,22 +33,22 @@ double& Material::Ci (const std::string& nameConst) { return MC[i]; } } - LOG(ERROR) << "Dan't find a material constant with name " << nameConst; - double dummy; + LOG(ERROR) << "Dan't find a material constant with name " << nameConst; + static double dummy; return dummy; } void Material::register_mat_const(uint16 num, ...) { - numC = num; - va_list vlist; - va_start(vlist, num); + numC = num; + va_list vlist; + va_start(vlist, num); MC_names.clear(); MC_names.reserve(numC); - for (uint16 i=0; i < num; i++) { - MC_names.push_back(va_arg(vlist,char*)); - } - MC = new double[num]; + for (uint16 i=0; i < num; i++) { + MC_names.push_back(va_arg(vlist,char*)); + } + MC = new double[num]; } } // namespace nla3d diff --git a/src/lib/materials/material.h b/src/lib/materials/material.h index c8b0d00..0b1cdb0 100755 --- a/src/lib/materials/material.h +++ b/src/lib/materials/material.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once #include @@ -14,57 +14,57 @@ class MaterialFactory; // Material must have named material constants class Material { public: - Material () : MC(NULL), numC(0), code(0) { - } - Material (uint16 num_c); - ~Material() // TODO: discover the virtual destructor - { - if (MC) delete[] MC; - MC = NULL; - } + Material () : MC(NULL), numC(0), code(0) { + } + Material (uint16 num_c); + ~Material() // TODO: discover the virtual destructor + { + if (MC) delete[] MC; + MC = NULL; + } - virtual std::string toString(); + virtual std::string toString(); - std::string getName(); - uint16 getCode (); + std::string getName(); + uint16 getCode (); // constants getters - double& Ci (uint16 i); - double& Ci (const std::string& nameConst); - uint16 getNumC (); + double& Ci (uint16 i); + double& Ci (const std::string& nameConst); + uint16 getNumC (); - static const double I[6]; + static const double I[6]; friend class MaterialFactory; protected: - void register_mat_const(uint16 num, ...); - double* MC; - std::vector MC_names; - uint16 numC; - uint16 code; - std::string name; + void register_mat_const(uint16 num, ...); + double* MC; + std::vector MC_names; + uint16 numC; + uint16 code; + std::string name; }; // ---=== FUNCTIONS ===--- // inline std::string Material::getName() { - return name; + return name; } inline uint16 Material::getCode () { - return code; + return code; } inline double& Material::Ci (uint16 i) { - assert(i < numC); - return MC[i]; + assert(i < numC); + return MC[i]; } inline uint16 Material::getNumC () { - return numC; + return numC; } } // namespace nla3d diff --git a/src/lib/materials/materials_hyperelastic.cpp b/src/lib/materials/materials_hyperelastic.cpp index 443ee8a..904e85b 100644 --- a/src/lib/materials/materials_hyperelastic.cpp +++ b/src/lib/materials/materials_hyperelastic.cpp @@ -1,30 +1,30 @@ // 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 +// https://github.com/dmitryikh/nla3d #include "materials/materials_hyperelastic.h" namespace nla3d { using namespace solidmech; -const double Mat_Hyper_Isotrop_General::II[6][6] = { {1.0,0.0,0.0,0.0,0.0,0.0}, - {0.0,0.5,0.0,0.0,0.0,0.0}, - {0.0,0.0,0.5,0.0,0.0,0.0}, - {0.0,0.0,0.0,1.0,0.0,0.0}, - {0.0,0.0,0.0,0.0,0.5,0.0}, - {0.0,0.0,0.0,0.0,0.0,1.0} }; +const double Mat_Hyper_Isotrop_General::II[6][6] = { {1.0,0.0,0.0,0.0,0.0,0.0}, + {0.0,0.5,0.0,0.0,0.0,0.0}, + {0.0,0.0,0.5,0.0,0.0,0.0}, + {0.0,0.0,0.0,1.0,0.0,0.0}, + {0.0,0.0,0.0,0.0,0.5,0.0}, + {0.0,0.0,0.0,0.0,0.0,1.0} }; //--------------------------------------------------------- //---------------Mat_Hyper_Isotrop_General----------------- //--------------------------------------------------------- void Mat_Hyper_Isotrop_General::getS_U (uint16 ncomp, const solidmech::tensorComponents* comps, const double* C, double *S) { - LOG(FATAL) << "Not now"; -} + LOG(FATAL) << "Not now"; +} void Mat_Hyper_Isotrop_General::getD_U (uint16 ncomp, const solidmech::tensorComponents* comps, const double* C, double *D) { - LOG(FATAL) << "Not now"; -} - + LOG(FATAL) << "Not now"; +} + //C - all times we must have 6 components void Mat_Hyper_Isotrop_General::getS_UP (uint16 ncomp, const solidmech::tensorComponents* comps, const double* C, const double press, double *S) { @@ -33,99 +33,99 @@ void Mat_Hyper_Isotrop_General::getS_UP (uint16 ncomp, const solidmech::tensorC double IC[3]; solidmech::IC_C(C, IC); - double _13I1C = 1.0/3.0*IC[0]; - double _23I2C = 2.0/3.0*IC[1]; - - double J = IC[2]; - - double oo = 1.0/(J*J); - double pp = pow(J,-2.0/3.0); - double pppp = pp*pp; - - double C_inv[6]; + double _13I1C = 1.0/3.0*IC[0]; + double _23I2C = 2.0/3.0*IC[1]; + + double J = IC[2]; + + double oo = 1.0/(J*J); + double pp = pow(J,-2.0/3.0); + double pppp = pp*pp; + + double C_inv[6]; solidmech::invC_C(C, J, C_inv); - - W_first_derivatives(IC[0]*pp, IC[1]*pppp, 0.0, alpha); //only first derivatives - double ko1 = 2*alpha[AL_1]*pp; - double ko2 = 2*alpha[AL_2]*pppp; - - double A[] = {I[0]-_13I1C*C_inv[0], I[1]-_13I1C*C_inv[1], I[2]-_13I1C*C_inv[2], - I[3]-_13I1C*C_inv[3], I[4]-_13I1C*C_inv[4], I[5]-_13I1C*C_inv[5]}; - double B[] = {IC[0]*I[0]-C[0]-_23I2C*C_inv[0], IC[0]*I[1]-C[1]-_23I2C*C_inv[1], IC[0]*I[2]-C[2]-_23I2C*C_inv[2], - IC[0]*I[3]-C[3]-_23I2C*C_inv[3], IC[0]*I[4]-C[4]-_23I2C*C_inv[4], IC[0]*I[5]-C[5]-_23I2C*C_inv[5]}; - solidmech::tensorComponents ij; - for (uint16 i=0; i < ncomp; i++) - { - ij = comps[i]; - S[i] = ko1*A[ij]+ko2*B[ij]+press*J*C_inv[ij]; - } + + W_first_derivatives(IC[0]*pp, IC[1]*pppp, 0.0, alpha); //only first derivatives + double ko1 = 2*alpha[AL_1]*pp; + double ko2 = 2*alpha[AL_2]*pppp; + + double A[] = {I[0]-_13I1C*C_inv[0], I[1]-_13I1C*C_inv[1], I[2]-_13I1C*C_inv[2], + I[3]-_13I1C*C_inv[3], I[4]-_13I1C*C_inv[4], I[5]-_13I1C*C_inv[5]}; + double B[] = {IC[0]*I[0]-C[0]-_23I2C*C_inv[0], IC[0]*I[1]-C[1]-_23I2C*C_inv[1], IC[0]*I[2]-C[2]-_23I2C*C_inv[2], + IC[0]*I[3]-C[3]-_23I2C*C_inv[3], IC[0]*I[4]-C[4]-_23I2C*C_inv[4], IC[0]*I[5]-C[5]-_23I2C*C_inv[5]}; + solidmech::tensorComponents ij; + for (uint16 i=0; i < ncomp; i++) + { + ij = comps[i]; + S[i] = ko1*A[ij]+ko2*B[ij]+press*J*C_inv[ij]; + } } void Mat_Hyper_Isotrop_General::getDdDp_UP (uint16 ncomp, const solidmech::tensorComponents* comps, const double* C, const double press, double *Dd, double *Dp) { - double alpha[5]; + double alpha[5]; double IC[3]; solidmech::IC_C(C, IC); - double _13I1C = 1.0/3.0*IC[0]; - double _23I2C = 2.0/3.0*IC[1]; - - double J = IC[2]; - - double pp = pow(J,-2.0/3.0); - double oo = pow(J,-2); - double pppp = pp*pp; - - double C_inv[6]; + double _13I1C = 1.0/3.0*IC[0]; + double _23I2C = 2.0/3.0*IC[1]; + + double J = IC[2]; + + double pp = pow(J,-2.0/3.0); + double oo = pow(J,-2); + double pppp = pp*pp; + + double C_inv[6]; solidmech::invC_C(C, J, C_inv); - /* - AL_1 = 0 - AL_2 = 1 - AL_11 = 2 - AL_12 = 3 - AL_22 = 4 - */ - W_second_derivatives(IC[0]*pp, IC[1]*pppp, 1.0, alpha); - - double A[] = {I[0]-_13I1C*C_inv[0], I[1]-_13I1C*C_inv[1], I[2]-_13I1C*C_inv[2], - I[3]-_13I1C*C_inv[3], I[4]-_13I1C*C_inv[4], I[5]-_13I1C*C_inv[5]}; - double B[] = {IC[0]*I[0]-C[0]-_23I2C*C_inv[0], IC[0]*I[1]-C[1]-_23I2C*C_inv[1], IC[0]*I[2]-C[2]-_23I2C*C_inv[2], - IC[0]*I[3]-C[3]-_23I2C*C_inv[3], IC[0]*I[4]-C[4]-_23I2C*C_inv[4], IC[0]*I[5]-C[5]-_23I2C*C_inv[5]}; - - double IIt[6][6] = {{-C_inv[M_XX]*C_inv[M_XX], -C_inv[M_XX]*C_inv[M_XY], -C_inv[M_XX]*C_inv[M_XZ], -C_inv[M_XY]*C_inv[M_XY], -C_inv[M_XY]*C_inv[M_XZ], -C_inv[M_XZ]*C_inv[M_XZ]}, - // 1211 = 11*12, 1212 = 0.5*(11*22 + 12*12), 1213 = 0.5*(11*23+13*12), 1222 = 12*22, 1223 = 0.5*(12*23+13*22), 1233 = 13*23 - {-C_inv[M_XX]*C_inv[M_XY], -0.5*(C_inv[M_XX]*C_inv[M_YY]+C_inv[M_XY]*C_inv[M_XY]), -0.5*(C_inv[M_XX]*C_inv[M_YZ]+C_inv[M_XZ]*C_inv[M_XY]), -C_inv[M_XY]*C_inv[M_YY], -0.5*(C_inv[M_XY]*C_inv[M_YZ]+C_inv[M_XZ]*C_inv[M_YY]), -C_inv[M_XZ]*C_inv[M_YZ]}, - // 1311 = 11*13, 1312 = 0.5*(11*23+12*13), 1313 = 0.5*(11*33+13*13), 1322 = 12*23, 1323 = 0.5*(12*33+13*23), 1333 = 13*33 - {-C_inv[M_XX]*C_inv[M_XZ], -0.5*(C_inv[M_XX]*C_inv[M_YZ]+C_inv[M_XY]*C_inv[M_XZ]), -0.5*(C_inv[M_XX]*C_inv[M_ZZ]+C_inv[M_XZ]*C_inv[M_XZ]), -C_inv[M_XY]*C_inv[M_YZ], -0.5*(C_inv[M_XY]*C_inv[M_ZZ]+C_inv[M_XZ]*C_inv[M_YZ]), -C_inv[M_XZ]*C_inv[M_ZZ]}, - // 2211 = 12*12, 2212 = 12*22, 2213 = 12*23, 2222 = 22*22, 2223 = 22*23, 2233 = 23*23 - {-C_inv[M_XY]*C_inv[M_XY], -C_inv[M_XY]*C_inv[M_YY], -C_inv[M_XY]*C_inv[M_YZ], -C_inv[M_YY]*C_inv[M_YY], -C_inv[M_YY]*C_inv[M_YZ], -C_inv[M_YZ]*C_inv[M_YZ]}, - // 2311 = 12*13, 2312 = 0.5*(12*23+22*13), 2313 = 0.5*(12*33+23*13), 2322 = 22*23, 2323 = 0.5*(22*33+23*23), 2333 = 23*33 - {-C_inv[M_XY]*C_inv[M_XZ], -0.5*(C_inv[M_XY]*C_inv[M_YZ]+C_inv[M_YY]*C_inv[M_XZ]), -0.5*(C_inv[M_XY]*C_inv[M_ZZ]+C_inv[M_YZ]*C_inv[M_XZ]), -C_inv[M_YY]*C_inv[M_YZ], -0.5*(C_inv[M_YY]*C_inv[M_ZZ]+C_inv[M_YZ]*C_inv[M_YZ]), -C_inv[M_YZ]*C_inv[M_ZZ]}, - // 3311 = 13*13, 3312 = 13*23, 3313 = 13*33, 3322 = 23*23, 3323 = 23*33, 3333 = 33*33 - {-C_inv[M_XZ]*C_inv[M_XZ], -C_inv[M_XZ]*C_inv[M_YZ], -C_inv[M_XZ]*C_inv[M_ZZ], -C_inv[M_YZ]*C_inv[M_YZ], -C_inv[M_YZ]*C_inv[M_ZZ], -C_inv[M_ZZ]*C_inv[M_ZZ]} - }; - - uint16 ind = 0; - solidmech::tensorComponents ij; - solidmech::tensorComponents kl; - for (uint16 i=0; i < ncomp; i++) - { - ij = comps[i]; - for (uint16 j=i; j < ncomp; j++) - { - - kl = comps[j]; - Dd[ind]=4.0*alpha[AL_11]*pppp*A[ij]*A[kl]+4.0*alpha[AL_12]*oo*(B[ij]*A[kl]+A[ij]*B[kl])+4.0*alpha[AL_22]*pppp*pppp*B[ij]*B[kl] - - 4.0/3.0*alpha[AL_1]*pp*(C_inv[ij]*A[kl]+A[ij]*C_inv[kl]+_13I1C*C_inv[ij]*C_inv[kl]+IC[0]*IIt[ij][kl]) - - 8.0/3.0*alpha[AL_2]*pppp*(C_inv[ij]*B[kl]+B[ij]*C_inv[kl]+_23I2C*C_inv[ij]*C_inv[kl]-3.0/2.0*I[ij]*I[kl]+3.0/2.0*II[ij][kl]+IC[1]*IIt[ij][kl])+ - press*J*(C_inv[ij]*C_inv[kl]+2*IIt[ij][kl]); - Dd[ind]=Dd[ind]/2.0; - ind++; - } - Dp[i] = J*C_inv[ij]; - } + /* + AL_1 = 0 + AL_2 = 1 + AL_11 = 2 + AL_12 = 3 + AL_22 = 4 + */ + W_second_derivatives(IC[0]*pp, IC[1]*pppp, 1.0, alpha); + + double A[] = {I[0]-_13I1C*C_inv[0], I[1]-_13I1C*C_inv[1], I[2]-_13I1C*C_inv[2], + I[3]-_13I1C*C_inv[3], I[4]-_13I1C*C_inv[4], I[5]-_13I1C*C_inv[5]}; + double B[] = {IC[0]*I[0]-C[0]-_23I2C*C_inv[0], IC[0]*I[1]-C[1]-_23I2C*C_inv[1], IC[0]*I[2]-C[2]-_23I2C*C_inv[2], + IC[0]*I[3]-C[3]-_23I2C*C_inv[3], IC[0]*I[4]-C[4]-_23I2C*C_inv[4], IC[0]*I[5]-C[5]-_23I2C*C_inv[5]}; + + double IIt[6][6] = {{-C_inv[M_XX]*C_inv[M_XX], -C_inv[M_XX]*C_inv[M_XY], -C_inv[M_XX]*C_inv[M_XZ], -C_inv[M_XY]*C_inv[M_XY], -C_inv[M_XY]*C_inv[M_XZ], -C_inv[M_XZ]*C_inv[M_XZ]}, + // 1211 = 11*12, 1212 = 0.5*(11*22 + 12*12), 1213 = 0.5*(11*23+13*12), 1222 = 12*22, 1223 = 0.5*(12*23+13*22), 1233 = 13*23 + {-C_inv[M_XX]*C_inv[M_XY], -0.5*(C_inv[M_XX]*C_inv[M_YY]+C_inv[M_XY]*C_inv[M_XY]), -0.5*(C_inv[M_XX]*C_inv[M_YZ]+C_inv[M_XZ]*C_inv[M_XY]), -C_inv[M_XY]*C_inv[M_YY], -0.5*(C_inv[M_XY]*C_inv[M_YZ]+C_inv[M_XZ]*C_inv[M_YY]), -C_inv[M_XZ]*C_inv[M_YZ]}, + // 1311 = 11*13, 1312 = 0.5*(11*23+12*13), 1313 = 0.5*(11*33+13*13), 1322 = 12*23, 1323 = 0.5*(12*33+13*23), 1333 = 13*33 + {-C_inv[M_XX]*C_inv[M_XZ], -0.5*(C_inv[M_XX]*C_inv[M_YZ]+C_inv[M_XY]*C_inv[M_XZ]), -0.5*(C_inv[M_XX]*C_inv[M_ZZ]+C_inv[M_XZ]*C_inv[M_XZ]), -C_inv[M_XY]*C_inv[M_YZ], -0.5*(C_inv[M_XY]*C_inv[M_ZZ]+C_inv[M_XZ]*C_inv[M_YZ]), -C_inv[M_XZ]*C_inv[M_ZZ]}, + // 2211 = 12*12, 2212 = 12*22, 2213 = 12*23, 2222 = 22*22, 2223 = 22*23, 2233 = 23*23 + {-C_inv[M_XY]*C_inv[M_XY], -C_inv[M_XY]*C_inv[M_YY], -C_inv[M_XY]*C_inv[M_YZ], -C_inv[M_YY]*C_inv[M_YY], -C_inv[M_YY]*C_inv[M_YZ], -C_inv[M_YZ]*C_inv[M_YZ]}, + // 2311 = 12*13, 2312 = 0.5*(12*23+22*13), 2313 = 0.5*(12*33+23*13), 2322 = 22*23, 2323 = 0.5*(22*33+23*23), 2333 = 23*33 + {-C_inv[M_XY]*C_inv[M_XZ], -0.5*(C_inv[M_XY]*C_inv[M_YZ]+C_inv[M_YY]*C_inv[M_XZ]), -0.5*(C_inv[M_XY]*C_inv[M_ZZ]+C_inv[M_YZ]*C_inv[M_XZ]), -C_inv[M_YY]*C_inv[M_YZ], -0.5*(C_inv[M_YY]*C_inv[M_ZZ]+C_inv[M_YZ]*C_inv[M_YZ]), -C_inv[M_YZ]*C_inv[M_ZZ]}, + // 3311 = 13*13, 3312 = 13*23, 3313 = 13*33, 3322 = 23*23, 3323 = 23*33, 3333 = 33*33 + {-C_inv[M_XZ]*C_inv[M_XZ], -C_inv[M_XZ]*C_inv[M_YZ], -C_inv[M_XZ]*C_inv[M_ZZ], -C_inv[M_YZ]*C_inv[M_YZ], -C_inv[M_YZ]*C_inv[M_ZZ], -C_inv[M_ZZ]*C_inv[M_ZZ]} + }; + + uint16 ind = 0; + solidmech::tensorComponents ij; + solidmech::tensorComponents kl; + for (uint16 i=0; i < ncomp; i++) + { + ij = comps[i]; + for (uint16 j=i; j < ncomp; j++) + { + + kl = comps[j]; + Dd[ind]=4.0*alpha[AL_11]*pppp*A[ij]*A[kl]+4.0*alpha[AL_12]*oo*(B[ij]*A[kl]+A[ij]*B[kl])+4.0*alpha[AL_22]*pppp*pppp*B[ij]*B[kl] - + 4.0/3.0*alpha[AL_1]*pp*(C_inv[ij]*A[kl]+A[ij]*C_inv[kl]+_13I1C*C_inv[ij]*C_inv[kl]+IC[0]*IIt[ij][kl]) - + 8.0/3.0*alpha[AL_2]*pppp*(C_inv[ij]*B[kl]+B[ij]*C_inv[kl]+_23I2C*C_inv[ij]*C_inv[kl]-3.0/2.0*I[ij]*I[kl]+3.0/2.0*II[ij][kl]+IC[1]*IIt[ij][kl])+ + press*J*(C_inv[ij]*C_inv[kl]+2*IIt[ij][kl]); + Dd[ind]=Dd[ind]/2.0; + ind++; + } + Dp[i] = J*C_inv[ij]; + } } @@ -134,17 +134,17 @@ void Mat_Hyper_Isotrop_General::getDdDp_UP (uint16 ncomp, const solidmech::tens //--------------------------------------------------------- void Mat_Comp_Neo_Hookean::W_first_derivatives (double I1, double I2, double I3, double* alpha) { - alpha[AL_1] = 0.5*MC[C_G]; - alpha[AL_2] = 0.0; + alpha[AL_1] = 0.5*MC[C_G]; + alpha[AL_2] = 0.0; } void Mat_Comp_Neo_Hookean::W_second_derivatives (double I1, double I2, double I3, double* alpha) { - alpha[AL_1] = 0.5*MC[C_G]; - alpha[AL_2] = 0.0; - alpha[AL_11] = 0.0; - alpha[AL_22] = 0.0; - alpha[AL_12] = 0.0; + alpha[AL_1] = 0.5*MC[C_G]; + alpha[AL_2] = 0.0; + alpha[AL_11] = 0.0; + alpha[AL_22] = 0.0; + alpha[AL_12] = 0.0; } double Mat_Comp_Neo_Hookean::W (double I1, double I2, double I3) { @@ -159,16 +159,16 @@ double Mat_Comp_Neo_Hookean::getK() { //-------------------Mat_Comp_Biderman--------------------- //--------------------------------------------------------- void Mat_Comp_Biderman::W_first_derivatives (double I1, double I2, double I3, double* alpha) { - alpha[AL_1] = MC[C_C10]+2*MC[C_C20]*(I1-3.0)+3*MC[C_C30]*(I1-3.0)*(I1-3.0); - alpha[AL_2] = MC[C_C01]; + alpha[AL_1] = MC[C_C10]+2*MC[C_C20]*(I1-3.0)+3*MC[C_C30]*(I1-3.0)*(I1-3.0); + alpha[AL_2] = MC[C_C01]; } void Mat_Comp_Biderman::W_second_derivatives (double I1, double I2, double I3, double* alpha) { - alpha[AL_1] = MC[C_C10]+2*MC[C_C20]*(I1-3.0)+3*MC[C_C30]*(I1-3.0)*(I1-3.0); - alpha[AL_2] = MC[C_C01]; - alpha[AL_11] = 2*MC[C_C20]+6*MC[C_C30]*(I1-3.0); - alpha[AL_22] = 0.0; - alpha[AL_12] = 0.0; + alpha[AL_1] = MC[C_C10]+2*MC[C_C20]*(I1-3.0)+3*MC[C_C30]*(I1-3.0)*(I1-3.0); + alpha[AL_2] = MC[C_C01]; + alpha[AL_11] = 2*MC[C_C20]+6*MC[C_C30]*(I1-3.0); + alpha[AL_22] = 0.0; + alpha[AL_12] = 0.0; } double Mat_Comp_Biderman::W (double I1, double I2, double I3) { return MC[C_C10]*(I1-3.0) + MC[C_C20]*(I1-3.0)*(I1-3.0) + MC[C_C30]*(I1-3.0)*(I1-3.0)*(I1-3.0) + MC[C_C01]*(I2-3.0); @@ -183,17 +183,17 @@ double Mat_Comp_Biderman::getK() { //--------------------------------------------------------- void Mat_Comp_MooneyRivlin::W_first_derivatives (double I1, double I2, double I3, double* alpha) { - alpha[AL_1] = MC[C_C10]; - alpha[AL_2] = MC[C_C01]; + alpha[AL_1] = MC[C_C10]; + alpha[AL_2] = MC[C_C01]; } void Mat_Comp_MooneyRivlin::W_second_derivatives (double I1, double I2, double I3, double* alpha) { - alpha[AL_1] = MC[C_C10]; - alpha[AL_2] = MC[C_C01]; - alpha[AL_11] = 0.0; - alpha[AL_22] = 0.0; - alpha[AL_12] = 0.0; + alpha[AL_1] = MC[C_C10]; + alpha[AL_2] = MC[C_C01]; + alpha[AL_11] = 0.0; + alpha[AL_22] = 0.0; + alpha[AL_12] = 0.0; } double Mat_Comp_MooneyRivlin::W (double I1, double I2, double I3) { return MC[C_C10]*(I1-3.0) + MC[C_C01]*(I2-3.0); diff --git a/src/lib/materials/materials_hyperelastic.h b/src/lib/materials/materials_hyperelastic.h index 6b81b28..39d7b08 100644 --- a/src/lib/materials/materials_hyperelastic.h +++ b/src/lib/materials/materials_hyperelastic.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once #include "materials/material.h" @@ -14,92 +14,92 @@ namespace nla3d { //----------------------------------------------// class Mat_Hyper_Isotrop_General : public Material { - public: - //Mat_Hyper_Isotrop_General () { - //} + public: + //Mat_Hyper_Isotrop_General () { + //} // constants for derives in double array - enum mat_func_deriv { - AL_1 = 0, - AL_2 = 1, - AL_11 = 2, - AL_12 = 3, - AL_22 = 4 - }; - //for non linear U elements - void getS_U (uint16 ncomp, const solidmech::tensorComponents* comps, const double* C, double *S); - void getD_U (uint16 ncomp, const solidmech::tensorComponents* comps, const double* C, double *D); - //for nonlinear U-P elements - void getS_UP (uint16 ncomp, const solidmech::tensorComponents* comps, const double* C, const double press, double *S); + enum mat_func_deriv { + AL_1 = 0, + AL_2 = 1, + AL_11 = 2, + AL_12 = 3, + AL_22 = 4 + }; + //for non linear U elements + void getS_U (uint16 ncomp, const solidmech::tensorComponents* comps, const double* C, double *S); + void getD_U (uint16 ncomp, const solidmech::tensorComponents* comps, const double* C, double *D); + //for nonlinear U-P elements + void getS_UP (uint16 ncomp, const solidmech::tensorComponents* comps, const double* C, const double press, double *S); void getDdDp_UP (uint16 ncomp, const solidmech::tensorComponents* comps, const double* C, const double press, double *Dd, double *Dp); - virtual void W_first_derivatives (double I1, double I2, double I3, double *alpha) = 0; - virtual void W_second_derivatives (double I1, double I2, double I3, double *alpha) = 0; + virtual void W_first_derivatives (double I1, double I2, double I3, double *alpha) = 0; + virtual void W_second_derivatives (double I1, double I2, double I3, double *alpha) = 0; virtual double W (double I1, double I2, double I3) = 0; virtual double getK() = 0; - - static const double II[6][6]; + + static const double II[6][6]; }; class Mat_Comp_Neo_Hookean : public Mat_Hyper_Isotrop_General { - public: - enum Mat_Comp_Neo_Hookean_CONST { - C_G = 0, - C_K - }; - - Mat_Comp_Neo_Hookean () { - register_mat_const(2,"G","K"); + public: + enum Mat_Comp_Neo_Hookean_CONST { + C_G = 0, + C_K + }; + + Mat_Comp_Neo_Hookean () { + register_mat_const(2,"G","K"); name = "Neo-Hookean"; - } + } - void W_first_derivatives (double I1, double I2, double I3, double *alpha); - void W_second_derivatives (double I1, double I2, double I3, double *alpha); - double W (double I1, double I2, double I3); + void W_first_derivatives (double I1, double I2, double I3, double *alpha); + void W_second_derivatives (double I1, double I2, double I3, double *alpha); + double W (double I1, double I2, double I3); double getK(); }; class Mat_Comp_Biderman : public Mat_Hyper_Isotrop_General { - public: - enum Mat_Comp_Biderman_CONST { - C_C10 = 0, - C_C20, - C_C30, - C_C01, - C_K - }; - - Mat_Comp_Biderman () { - register_mat_const(5,"C10","C20","C30","C01","K"); + public: + enum Mat_Comp_Biderman_CONST { + C_C10 = 0, + C_C20, + C_C30, + C_C01, + C_K + }; + + Mat_Comp_Biderman () { + register_mat_const(5,"C10","C20","C30","C01","K"); name = "Biderman"; - } - void W_first_derivatives (double I1, double I2, double I3, double *alpha); - void W_second_derivatives (double I1, double I2, double I3, double *alpha); - double W (double I1, double I2, double I3); + } + void W_first_derivatives (double I1, double I2, double I3, double *alpha); + void W_second_derivatives (double I1, double I2, double I3, double *alpha); + double W (double I1, double I2, double I3); double getK(); }; class Mat_Comp_MooneyRivlin : public Mat_Hyper_Isotrop_General { - public: - enum Mat_Comp_MooneyRivlin_CONST { - C_C10 = 0, - C_C01, - C_K - }; - - Mat_Comp_MooneyRivlin () { - register_mat_const(3,"C10","C01","K"); + public: + enum Mat_Comp_MooneyRivlin_CONST { + C_C10 = 0, + C_C01, + C_K + }; + + Mat_Comp_MooneyRivlin () { + register_mat_const(3,"C10","C01","K"); name = "Money-Rivlin"; - } - void W_first_derivatives (double I1, double I2, double I3, double *alpha); - void W_second_derivatives (double I1, double I2, double I3, double *alpha); - double W (double I1, double I2, double I3); + } + void W_first_derivatives (double I1, double I2, double I3, double *alpha); + void W_second_derivatives (double I1, double I2, double I3, double *alpha); + double W (double I1, double I2, double I3); double getK(); }; diff --git a/src/lib/math/BlockSparseMatrix.h b/src/lib/math/BlockSparseMatrix.h index b524975..b9668f6 100644 --- a/src/lib/math/BlockSparseMatrix.h +++ b/src/lib/math/BlockSparseMatrix.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once @@ -15,7 +15,7 @@ namespace math { // global sparse matrix is represented by different blocks (`nb` is number of blocks). Here is a // illustration: // Here is global symmetric matrix A `n` by `n`, it consists of 3 distinct sparse matrix blocks: -// | A1 | A12 | +// | A1 | A12 | // A = |-----------|, // |A12^T| A2 | // @@ -41,8 +41,8 @@ class BlockSparseSymMatrix { public: BlockSparseSymMatrix(std::initializer_list _rows_in_block, uint32 max_in_row = 100); BlockSparseSymMatrix(BlockSparseSymMatrix* ex); - - // add non-zero entry to sparse matrix. This should be called before compress(). + + // add non-zero entry to sparse matrix. This should be called before compress(). void addEntry(uint32 _i, uint32 _j); // add value to the _i, _j entry. This should be called after compress(). diff --git a/src/lib/math/CMakeLists.txt b/src/lib/math/CMakeLists.txt index f171312..73f2b0b 100755 --- a/src/lib/math/CMakeLists.txt +++ b/src/lib/math/CMakeLists.txt @@ -1,4 +1,4 @@ -# Collect sources into the variable MATH_SOURCES with +# Collect sources into the variable MATH_SOURCES with # having to explicitly list each header and source file file (GLOB MATH_SOURCES "*.cpp") diff --git a/src/lib/math/EquationSolver.cpp b/src/lib/math/EquationSolver.cpp index d0e8e43..ad7786d 100644 --- a/src/lib/math/EquationSolver.cpp +++ b/src/lib/math/EquationSolver.cpp @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #include "math/EquationSolver.h" @@ -9,7 +9,6 @@ #endif //NLA3D_USE_MKL namespace nla3d { - namespace math { EquationSolver* defaultEquationSolver = new GaussDenseEquationSolver; @@ -43,7 +42,7 @@ void GaussDenseEquationSolver::substituteEquations(math::SparseSymMatrix* matrix // because this is simplest solver dedicated to perform functional tests // CHECK(nrhs == 1) << "GaussDenseEquationSolver support only 1 set of rhs values"; - + if (matA.dM() == nEq && matA.dN() == nEq) { matA.zero(); @@ -149,18 +148,18 @@ void PARDISO_equationSolver::solveEquations (math::SparseSymMatrix* matrix, doub void PARDISO_equationSolver::substituteEquations(math::SparseSymMatrix* matrix, double* rhs, double* unknowns) { - //Back substitution and iterative refinement + //Back substitution and iterative refinement // initialize error code - int error = 0; - int phase = 33; + int error = 0; + int phase = 33; int n = static_cast (nEq); - PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, matrix->getValuesArray(), + PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, matrix->getValuesArray(), (int*) matrix->getIofeirArray(), (int*) matrix->getColumnsArray(), - NULL, &nrhs, iparm, &msglvl, rhs, unknowns, &error); + NULL, &nrhs, iparm, &msglvl, rhs, unknowns, &error); - CHECK(error == 0) << "ERROR during solution. Error code = " << error; + CHECK(error == 0) << "ERROR during solution. Error code = " << error; } @@ -172,26 +171,26 @@ void PARDISO_equationSolver::factorizeEquations(math::SparseSymMatrix* matrix) { firstRun = false; CHECK(nEq == matrix->nRows()); - - int phase; + + int phase; // initialize error code - int error = 0; + int error = 0; - // phase 22 is the numerical factorization - phase = 22; + // phase 22 is the numerical factorization + phase = 22; int n = static_cast (nEq); - PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, matrix->getValuesArray(), + PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, matrix->getValuesArray(), (int*) matrix->getIofeirArray(), (int*) matrix->getColumnsArray(), - NULL, &nrhs, iparm, &msglvl, NULL, NULL, &error); + NULL, &nrhs, iparm, &msglvl, NULL, NULL, &error); CHECK(error == 0) << "ERROR during numerical factorization. Error code = " << error; } void PARDISO_equationSolver::initializePARDISO (math::SparseSymMatrix* matrix) { - for (uint16 i = 0; i < 64; i++) { + for (uint16 i = 0; i < 64; i++) { iparm[i]=0; } @@ -199,23 +198,23 @@ void PARDISO_equationSolver::initializePARDISO (math::SparseSymMatrix* matrix) { pt[i]=0; } - iparm[0] = 1; //no solver default - iparm[1] = 2; //fill-in reordering from meris - iparm[2] = MKL_Get_Max_Threads(); - iparm[3] = 0; //no iterative-direct algorithm - iparm[4] = 0; //no user fill-in reducing permutation - iparm[5] = 0; //write solution into x - iparm[6] = 16; //default logical fortran unit number for output - iparm[7] = 2; //max numbers of iterative refinement steps - iparm[9] = 13; //pertrub the pivor elements with 1e-13 - iparm[10] = 1; //use nonsymmetric permutation and scaling MPS - iparm[13]=0; //output: number of perturbed pivots - iparm[17]=-1; //output: number of nonzeros in the factor LU - iparm[18]=-1; //output: MFLOPS for LU factorization - iparm[19] = 0; //output: number of CG Iterations + iparm[0] = 1; //no solver default + iparm[1] = 2; //fill-in reordering from meris + iparm[2] = MKL_Get_Max_Threads(); + iparm[3] = 0; //no iterative-direct algorithm + iparm[4] = 0; //no user fill-in reducing permutation + iparm[5] = 0; //write solution into x + iparm[6] = 16; //default logical fortran unit number for output + iparm[7] = 2; //max numbers of iterative refinement steps + iparm[9] = 13; //pertrub the pivor elements with 1e-13 + iparm[10] = 1; //use nonsymmetric permutation and scaling MPS + iparm[13]=0; //output: number of perturbed pivots + iparm[17]=-1; //output: number of nonzeros in the factor LU + iparm[18]=-1; //output: MFLOPS for LU factorization + iparm[19] = 0; //output: number of CG Iterations LOG_IF(!isSymmetric, FATAL) << "For now PARDISO_equationSolver doesn't support non-symmetric matrices"; - if (isPositive) { + if (isPositive) { mtype = 2; LOG(INFO) << "EquationSolver will use positive symmetric solver"; } else { @@ -227,14 +226,14 @@ void PARDISO_equationSolver::initializePARDISO (math::SparseSymMatrix* matrix) { int n = static_cast (nEq); // initialize error code - int error = 0; + int error = 0; int phase = 11; PARDISO(pt, &maxfct, &mnum, &mtype,&phase, &n, matrix->getValuesArray(), - (int*) matrix->getIofeirArray(), - (int*) matrix->getColumnsArray(), - NULL, &nrhs, iparm, &msglvl, NULL, NULL, &error); + (int*) matrix->getIofeirArray(), + (int*) matrix->getColumnsArray(), + NULL, &nrhs, iparm, &msglvl, NULL, NULL, &error); CHECK(error == 0) << "ERROR during symbolic factorization. Error code = " << error; LOG(INFO) << "Number of nonzeros in factors = " << iparm[17] << ", number of factorization MFLOPS = " << iparm[18]; } @@ -244,13 +243,12 @@ void PARDISO_equationSolver::releasePARDISO () { int n = static_cast (nEq); // initialize error code - int error = 0; - //Termination and release of memory - PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, NULL, NULL, NULL, NULL, &nrhs, iparm, &msglvl, NULL, NULL, &error); + int error = 0; + //Termination and release of memory + PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, NULL, NULL, NULL, NULL, &nrhs, iparm, &msglvl, NULL, NULL, &error); LOG_IF (error != 0, WARNING) << "ERROR during PARDISO termination. Error code = " << error; } #endif //NLA3D_USE_MKL } //namespace math - } //namespace nla3d diff --git a/src/lib/math/EquationSolver.h b/src/lib/math/EquationSolver.h index f0f3e89..d6dd5ef 100644 --- a/src/lib/math/EquationSolver.h +++ b/src/lib/math/EquationSolver.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once #include "sys.h" @@ -8,7 +8,6 @@ #include "math/SparseMatrix.h" namespace nla3d { - namespace math { class SparseSymMatrix; @@ -27,8 +26,8 @@ class EquationSolver { protected: uint32 nEq = 0; - // number of rhs - int nrhs = 1; + // number of rhs + int nrhs = 1; bool isSymmetric = true; bool isPositive = true; }; @@ -57,26 +56,25 @@ class PARDISO_equationSolver : public EquationSolver { void releasePARDISO (); // Internal solver memory pointer pt - void *pt[64]; + void *pt[64]; // Paramaters for PARDISO solver (see MKL manual for clarifications) - int iparm[64]; + int iparm[64]; // maximum number of numerical factorizations - int maxfct = 1; + int maxfct = 1; // which factorization to use - int mnum = 1; + int mnum = 1; // don't print statistical information in file - int msglvl = 0; - + int msglvl = 0; + bool firstRun = true; // real symmetric undifinite defined matrix - int mtype = -2; + int mtype = -2; }; #endif //NLA3D_USE_MKL extern EquationSolver* defaultEquationSolver; } // namespace math - } //namespace nla3d diff --git a/src/lib/math/Mat.cpp b/src/lib/math/Mat.cpp index faf0c46..170bd5e 100755 --- a/src/lib/math/Mat.cpp +++ b/src/lib/math/Mat.cpp @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #include "Mat.h" @@ -12,17 +12,17 @@ namespace math { template<> double Mat<1,1>::det() { - return data[0][0]; + return data[0][0]; } template<> double Mat<2,2>::det() { - return data[0][0]*data[1][1]-data[0][1]*data[1][0]; + return data[0][0]*data[1][1]-data[0][1]*data[1][0]; } template<> double Mat<3,3>::det() { - return data[0][0]*(data[1][1]*data[2][2]-data[1][2]*data[2][1])-data[0][1]*(data[1][0]*data[2][2]-data[1][2]*data[2][0])+data[0][2]*(data[1][0]*data[2][1]-data[1][1]*data[2][0]); + return data[0][0]*(data[1][1]*data[2][2]-data[1][2]*data[2][1])-data[0][1]*(data[1][0]*data[2][2]-data[1][2]*data[2][0])+data[0][2]*(data[1][0]*data[2][1]-data[1][1]*data[2][0]); } template<> @@ -33,29 +33,29 @@ Mat<1,1> Mat<1,1>::inv(double det) { template<> Mat<2,2> Mat<2,2>::inv(double det) { assert(det); - Mat<2,2> tmp(data[1][1],-data[0][1],-data[1][0],data[0][0]); - return tmp*(1.0/det); + Mat<2,2> tmp(data[1][1],-data[0][1],-data[1][0],data[0][0]); + return tmp*(1.0/det); } template<> Mat<3,3> Mat<3,3>::inv(double det) { assert(det); - Mat<3,3> tmp(data[1][1]*data[2][2]-data[1][2]*data[2][1],data[0][2]*data[2][1]-data[0][1]*data[2][2],data[0][1]*data[1][2]-data[0][2]*data[1][1], - data[2][0]*data[1][2]-data[1][0]*data[2][2],data[0][0]*data[2][2]-data[0][2]*data[2][0],data[1][0]*data[0][2]-data[0][0]*data[1][2], - data[1][0]*data[2][1]-data[1][1]*data[2][0],data[0][1]*data[2][0]-data[0][0]*data[2][1],data[0][0]*data[1][1]-data[0][1]*data[1][0]); - return tmp*(1.0/det); + Mat<3,3> tmp(data[1][1]*data[2][2]-data[1][2]*data[2][1],data[0][2]*data[2][1]-data[0][1]*data[2][2],data[0][1]*data[1][2]-data[0][2]*data[1][1], + data[2][0]*data[1][2]-data[1][0]*data[2][2],data[0][0]*data[2][2]-data[0][2]*data[2][0],data[1][0]*data[0][2]-data[0][0]*data[1][2], + data[1][0]*data[2][1]-data[1][1]*data[2][0],data[0][1]*data[2][0]-data[0][0]*data[2][1],data[0][0]*data[1][1]-data[0][1]*data[1][0]); + return tmp*(1.0/det); } //---------operator<<---------------------------------------------------------- std::ostream &operator<<(std::ostream &stream, dMat &obj) { - for (uint16 i = 0; i < obj.dimM; i++) - { - std::cout << "["; - for (uint16 j = 0; j < obj.dimN; j++) - std::cout << obj[i][j] << " "; - std::cout << "]" << endl; - } - return stream; + for (uint16 i = 0; i < obj.dimM; i++) + { + std::cout << "["; + for (uint16 j = 0; j < obj.dimN; j++) + std::cout << obj[i][j] << " "; + std::cout << "]" << endl; + } + return stream; } } // namespace math diff --git a/src/lib/math/Mat.h b/src/lib/math/Mat.h index cbb128e..67b6c48 100755 --- a/src/lib/math/Mat.h +++ b/src/lib/math/Mat.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once #include "sys.h" @@ -18,221 +18,221 @@ template class Mat { public: - Mat() { - assert(dimM && dimN); - } - Mat(double first, ...) { - assert(dimM && dimN); - va_list argp; - va_start(argp, first); - for (uint16 i=0; i < dimM; i++) - for (uint16 j=0; j < dimN; j++) - { - if (i == 0 && j == 0) data[i][j]=first; - else data[i][j]=va_arg(argp, double); - } - va_end(argp); - } - ~Mat() { } - - Vec& operator[] (uint16 n); - const Vec& operator[] (uint16 n) const; - Mat transpose (); - void display (); - Mat operator+ (const Mat &op); - Mat& operator+= (const Mat &op); - Mat operator- (); - Mat& operator= (const Mat &op); - Mat& operator= (const Vec &op); //TODO: CHECK this operation - Mat operator- (const Mat &op); - Mat operator* (const double op); - void Identity (); - void zero (); - double det(); - Vec eigenvalues(); - Mat inv(double det); - Mat cross_cut (uint16 cuti, uint16 cutj); - string toString (); - double* ptr (); - bool compare (const Mat &B, double eps = 0.00001); - void simple_read (std::istream &st); - //friend функции - template - friend std::ostream& operator<< (std::ostream& stream, const Mat &obj); - template + Mat() { + assert(dimM && dimN); + } + Mat(double first, ...) { + assert(dimM && dimN); + va_list argp; + va_start(argp, first); + for (uint16 i=0; i < dimM; i++) + for (uint16 j=0; j < dimN; j++) + { + if (i == 0 && j == 0) data[i][j]=first; + else data[i][j]=va_arg(argp, double); + } + va_end(argp); + } + ~Mat() { } + + Vec& operator[] (uint16 n); + const Vec& operator[] (uint16 n) const; + Mat transpose (); + void display (); + Mat operator+ (const Mat &op); + Mat& operator+= (const Mat &op); + Mat operator- (); + Mat& operator= (const Mat &op); + Mat& operator= (const Vec &op); //TODO: CHECK this operation + Mat operator- (const Mat &op); + Mat operator* (const double op); + void Identity (); + void zero (); + double det(); + Vec eigenvalues(); + Mat inv(double det); + Mat cross_cut (uint16 cuti, uint16 cutj); + string toString (); + double* ptr (); + bool compare (const Mat &B, double eps = 0.00001); + void simple_read (std::istream &st); + //friend функции + template + friend std::ostream& operator<< (std::ostream& stream, const Mat &obj); + template friend Mat operator* (const Mat &op1, const Mat &op2); - template + template friend Vec operator* (const Mat &op1, const Vec &op2); // template // friend bool matCompare (const Mat& mat1, const Mat& mat2, const double eps); - Vec data[dimM]; + Vec data[dimM]; private: //data was moved from here to public }; //---------operator<<---------------------------------------------------------- template std::ostream &operator<<(std::ostream &stream, const Mat &obj) { - for (uint16 i = 0; i < dimM; i++) { - stream << "[" << obj.data[i] << "]" << std::endl; - } - return stream; + for (uint16 i = 0; i < dimM; i++) { + stream << "[" << obj.data[i] << "]" << std::endl; + } + return stream; } //----------operator[]--------------------------------------------------------- template inline Vec& Mat::operator[] (uint16 n) { - assert(n < dimM); - return data[n]; + assert(n < dimM); + return data[n]; } template inline const Vec& Mat::operator[] (uint16 n) const { - assert(n < dimM); - return data[n]; + assert(n < dimM); + return data[n]; } //----------Identity ()--------------------------------------------------------- template void Mat::Identity () { - assert (dimM==dimN); - for (uint16 i=0; i < dimM; i++) - for (uint16 j=0; j < dimN; j++) - if (i==j) - data[i][j]=1.0f; - else - data[i][j]=0.0f; + assert (dimM==dimN); + for (uint16 i=0; i < dimM; i++) + for (uint16 j=0; j < dimN; j++) + if (i==j) + data[i][j]=1.0f; + else + data[i][j]=0.0f; } //----------zero()--------------------------------------------------------- template void Mat::zero() { - memset(data,0,sizeof(double)*dimM*dimN); + memset(data,0,sizeof(double)*dimM*dimN); } //-------------display()------------------------------------------------------ template void Mat::display() { - for (unsigned int i = 0; i < dimM; i++) { - std::cout << "["; - data[i].display(); - std::cout << "]" << std::endl; - } + for (unsigned int i = 0; i < dimM; i++) { + std::cout << "["; + data[i].display(); + std::cout << "]" << std::endl; + } } //----------transpose()------------------------------------------------------- template Mat Mat::transpose () { - Mat p; - for (uint16 i=0; i < dimM; i++) - for (uint16 j=0; j < dimN; j++) - p[j][i]=this->data[i][j]; - return p; + Mat p; + for (uint16 i=0; i < dimM; i++) + for (uint16 j=0; j < dimN; j++) + p[j][i]=this->data[i][j]; + return p; } //----------operator+(Mat)--------------------------------------------------------- template Mat Mat::operator+ (const Mat &op) { - Mat p; - for (uint16 i=0; i < dimM; i++) - for (uint16 j=0; j < dimN; j++) - p[i][j]=this->data[i][j]+op.data[i][j]; - return p; + Mat p; + for (uint16 i=0; i < dimM; i++) + for (uint16 j=0; j < dimN; j++) + p[i][j]=this->data[i][j]+op.data[i][j]; + return p; } //----------operator-(Mat)--------------------------------------------------------- //template //Mat Mat::operator- (const Mat &op) { -// Mat p; -// for (uint16 i=0; i < dimM; i++) -// for (uint16 j=0; j < dimN; j++) -// p[i][j]=this->data[i][j]-op.data[i][j]; -// return p; +// Mat p; +// for (uint16 i=0; i < dimM; i++) +// for (uint16 j=0; j < dimN; j++) +// p[i][j]=this->data[i][j]-op.data[i][j]; +// return p; //} //-----------operator*(double)---------------------------------------------------------- template Mat Mat::operator*(const double op) { - Mat p; - for (uint16 i=0; i < dimM; i++) - for (uint16 j=0; j < dimN; j++) - p[i][j]=this->data[i][j]*op; - return p; + Mat p; + for (uint16 i=0; i < dimM; i++) + for (uint16 j=0; j < dimN; j++) + p[i][j]=this->data[i][j]*op; + return p; } //----------operator+=--------------------------------------------------------- template Mat& Mat::operator+= (const Mat &op) { - for (uint16 i=0; i < dimM; i++) - for (uint16 j=0; j < dimN; j++) - this->data[i][j]+=op.data[i][j]; - return *this; + for (uint16 i=0; i < dimM; i++) + for (uint16 j=0; j < dimN; j++) + this->data[i][j]+=op.data[i][j]; + return *this; } //----------operator=(Mat)--------------------------------------------------------- template Mat& Mat::operator= (const Mat &op) { - for (uint16 i = 0; i < dimM; i++) - this->data[i]=op.data[i]; - return *this; + for (uint16 i = 0; i < dimM; i++) + this->data[i]=op.data[i]; + return *this; } //----------operator=(Vec)--------------------------------------------------------- template Mat& Mat::operator= (const Vec &op) { - assert(dimN == 1); //только для матрицы-столбца - for (uint16 i = 0; i < dimM; i++) - this->data[i][0] = op[i]; - return *this; + assert(dimN == 1); //только для матрицы-столбца + for (uint16 i = 0; i < dimM; i++) + this->data[i][0] = op[i]; + return *this; } //----------operator*(Mat)--------------------------------------------------------- template Mat operator*(const Mat &op1, const Mat &op2) { - assert(dimN1 == dimM2); - Mat p; - double element; - double *ptr1 = (double*) op1.data; - double *ptr2 = (double*) op2.data; - for (uint16 i=0; i < dimM1; i++) - for (uint16 j=0; j < dimN2; j++) - { - element=0.0f; - for (uint16 l=0; l < dimN1; l++) element+=ptr1[i*dimN1+l]*ptr2[l*dimN2+j]; - //(op1.data[i])[l]*(op2.data[l])[j]; - p[i][j]=element; - } - return p; - - - //assert(dimN == dimM2); - //Mat p; - //double element; - //for (uint16 i=0; i < dimM; i++) - // for (uint16 j=0; j < dimN2; j++) - // { - // element=0.0f; - // for (uint16 l=0; l < dimN; l++) element+=(op1.data[i])[l]*(op2.data[l])[j]; - // p[i][j]=element; - // } - //return p; + assert(dimN1 == dimM2); + Mat p; + double element; + double *ptr1 = (double*) op1.data; + double *ptr2 = (double*) op2.data; + for (uint16 i=0; i < dimM1; i++) + for (uint16 j=0; j < dimN2; j++) + { + element=0.0f; + for (uint16 l=0; l < dimN1; l++) element+=ptr1[i*dimN1+l]*ptr2[l*dimN2+j]; + //(op1.data[i])[l]*(op2.data[l])[j]; + p[i][j]=element; + } + return p; + + + //assert(dimN == dimM2); + //Mat p; + //double element; + //for (uint16 i=0; i < dimM; i++) + // for (uint16 j=0; j < dimN2; j++) + // { + // element=0.0f; + // for (uint16 l=0; l < dimN; l++) element+=(op1.data[i])[l]*(op2.data[l])[j]; + // p[i][j]=element; + // } + //return p; } //-----------operator*(Vec)--------------------------------------------------------- template Vec operator* (const Mat &op1, const Vec &op2) { - assert(dimN1==dimM2); - Vec p; - double el; - for (uint16 i=0; i < dimM1; i++) - { - el = 0.0f; - for (uint16 j=0; j < dimN1; j++) - el += op1.data[i][j]*op2[j]; - p[i] = el; - } - return p; + assert(dimN1==dimM2); + Vec p; + double el; + for (uint16 i=0; i < dimM1; i++) + { + el = 0.0f; + for (uint16 j=0; j < dimN1; j++) + el += op1.data[i][j]*op2[j]; + p[i] = el; + } + return p; } //-------------------------------------------------------------------------------- template string Mat::toString() { - string p; - for (unsigned int i = 0; i < dimM; i++) { - p+= "["; - p+=data[i].toString(); - p+= "]"; - if (i!=dimM-1) p+= "\n"; - } - return p; + string p; + for (unsigned int i = 0; i < dimM; i++) { + p+= "["; + p+=data[i].toString(); + p+= "]"; + if (i!=dimM-1) p+= "\n"; + } + return p; } //------------------------------------------------------------------------------- template<> double Mat<1,1>::det(); @@ -242,17 +242,17 @@ template<> double Mat<3,3>::det(); template double Mat::det() { - //разложение по строке i=0 - assert(dimM == dimN); - double det = 0.0; - double koef = 1; - for (uint16 j=0; j < dimN; j++) - { - Mat tmp; - det += koef*data[0][j]*cross_cut(0,j).det(); - koef *= -1; - } - return det; + //разложение по строке i=0 + assert(dimM == dimN); + double det = 0.0; + double koef = 1; + for (uint16 j=0; j < dimN; j++) + { + Mat tmp; + det += koef*data[0][j]*cross_cut(0,j).det(); + koef *= -1; + } + return det; } template<> Mat<1, 1> Mat<1,1>::inv(double det); @@ -262,138 +262,138 @@ template<> Mat<3, 3> Mat<3,3>::inv(double det); template Mat Mat::inv(double det) { - //через матрицу алг. дополнений - assert(dimM == dimN); - Mat C; - for (uint16 i=0; i < dimM; i++) - for (uint16 j=0; j < dimN; j++) - C[i][j] = npow(-1,i+j)*cross_cut(i,j).det(); //CHECK - C = C.transpose()*(1/det); - return C; + //через матрицу алг. дополнений + assert(dimM == dimN); + Mat C; + for (uint16 i=0; i < dimM; i++) + for (uint16 j=0; j < dimN; j++) + C[i][j] = npow(-1,i+j)*cross_cut(i,j).det(); //CHECK + C = C.transpose()*(1/det); + return C; } template Mat Mat::cross_cut (uint16 cuti, uint16 cutj) { - Mat tmp; - for (uint16 ii=0; ii < dimM-1; ii++) - for (uint16 jj=0; jj < dimN-1; jj++) - tmp[ii][jj] = data[(ii>=cuti)?ii+1:ii][(jj>=cutj)?jj+1:jj]; - return tmp; + Mat tmp; + for (uint16 ii=0; ii < dimM-1; ii++) + for (uint16 jj=0; jj < dimN-1; jj++) + tmp[ii][jj] = data[(ii>=cuti)?ii+1:ii][(jj>=cutj)?jj+1:jj]; + return tmp; } template Vec Mat::eigenvalues() { - assert(dimM == 3 && dimN == 3); //TODO: пока только для матриц 3 на 3 - double I1 = data[0][0] + data[1][1] + data[2][2]; - double I2 = data[0][0]*data[1][1] + data[1][1]*data[2][2] + data[0][0]*data[2][2] - - data[0][1]*data[0][1] - data[0][2]*data[0][2] - data[1][2]*data[1][2]; - double I3 = data[0][0]*data[1][1]*data[2][2] + 2*data[0][1]*data[1][2]*data[0][2] - - data[0][0]*data[1][2]*data[1][2]-data[1][1]*data[0][2]*data[0][2]-data[2][2]*data[0][1]*data[0][1]; - - double a = -I1; - double b = I2; - double c = -I3; - // код из http://ru.wikipedia.org/wiki/%D0%A2%D1%80%D0%B8%D0%B3%D0%BE%D0%BD%D0%BE%D0%BC%D0%B5%D1%82%D1%80%D0%B8%D1%87%D0%B5%D1%81%D0%BA%D0%B0%D1%8F_%D1%84%D0%BE%D1%80%D0%BC%D1%83%D0%BB%D0%B0_%D0%92%D0%B8%D0%B5%D1%82%D0%B0 - // для трех вещественных корней - // x*x*x + a * x * x + b * x + c == 0 - double p = b - a * a / 3; - double q = 2 * a * a * a / 27 - a * b / 3 + c; - double A = sqrt(- 4 * p / 3); - - double c3phi = - 4 * q / (A * A * A); - - double phi = acos(c3phi) / 3; - - double root1 = A * cos(phi) - a / 3; - double root2 = A * cos(phi + 2 * M_PI / 3) - a / 3; - double root3 = A * cos(phi - 2 * M_PI / 3) - a / 3; - - //сортируем - double s1,s2,s3; - if (root1 > root2 && root1 > root3) - { - s1 = root1; - if (root2 > root3) - { - s2 = root2; - s3 = root3; - } - else - { - s2 = root3; - s3 = root2; - } - } - else if (root2 > root1 && root2 > root3) - { - s1 = root2; - if (root1 > root3) - { - s2 = root1; - s3 = root3; - } - else - { - s2 = root3; - s3 = root1; - } - } - else - { - s1 = root3; - if (root1 > root2) - { - s2 = root1; - s3 = root2; - } - else - { - s2 = root2; - s3 = root1; - } - } - return Vec<3>(s1,s2,s3); + assert(dimM == 3 && dimN == 3); //TODO: пока только для матриц 3 на 3 + double I1 = data[0][0] + data[1][1] + data[2][2]; + double I2 = data[0][0]*data[1][1] + data[1][1]*data[2][2] + data[0][0]*data[2][2] - + data[0][1]*data[0][1] - data[0][2]*data[0][2] - data[1][2]*data[1][2]; + double I3 = data[0][0]*data[1][1]*data[2][2] + 2*data[0][1]*data[1][2]*data[0][2] - + data[0][0]*data[1][2]*data[1][2]-data[1][1]*data[0][2]*data[0][2]-data[2][2]*data[0][1]*data[0][1]; + + double a = -I1; + double b = I2; + double c = -I3; + // код из http://ru.wikipedia.org/wiki/%D0%A2%D1%80%D0%B8%D0%B3%D0%BE%D0%BD%D0%BE%D0%BC%D0%B5%D1%82%D1%80%D0%B8%D1%87%D0%B5%D1%81%D0%BA%D0%B0%D1%8F_%D1%84%D0%BE%D1%80%D0%BC%D1%83%D0%BB%D0%B0_%D0%92%D0%B8%D0%B5%D1%82%D0%B0 + // для трех вещественных корней + // x*x*x + a * x * x + b * x + c == 0 + double p = b - a * a / 3; + double q = 2 * a * a * a / 27 - a * b / 3 + c; + double A = sqrt(- 4 * p / 3); + + double c3phi = - 4 * q / (A * A * A); + + double phi = acos(c3phi) / 3; + + double root1 = A * cos(phi) - a / 3; + double root2 = A * cos(phi + 2 * M_PI / 3) - a / 3; + double root3 = A * cos(phi - 2 * M_PI / 3) - a / 3; + + //сортируем + double s1,s2,s3; + if (root1 > root2 && root1 > root3) + { + s1 = root1; + if (root2 > root3) + { + s2 = root2; + s3 = root3; + } + else + { + s2 = root3; + s3 = root2; + } + } + else if (root2 > root1 && root2 > root3) + { + s1 = root2; + if (root1 > root3) + { + s2 = root1; + s3 = root3; + } + else + { + s2 = root3; + s3 = root1; + } + } + else + { + s1 = root3; + if (root1 > root2) + { + s2 = root1; + s3 = root2; + } + else + { + s2 = root2; + s3 = root1; + } + } + return Vec<3>(s1,s2,s3); } template double* Mat::ptr () { - return data[0].ptr(); + return data[0].ptr(); } template bool Mat::compare (const Mat &B, double eps) { - double *Dp = (double*) data; - // double *Bp = B.ptr(); + double *Dp = (double*) data; + // double *Bp = B.ptr(); double *Bp = (double*) B.data; - for (uint16 i = 0; i < dimM; i++) { - for (uint16 j = 0; j < dimN; j++) { - if (fabs(*Dp-*Bp) > eps) { + for (uint16 i = 0; i < dimM; i++) { + for (uint16 j = 0; j < dimN; j++) { + if (fabs(*Dp-*Bp) > eps) { LOG(INFO) << "Mat[" << i << "][" << j << "]: " << *Dp << " != " << *Bp; - return false; - } - Dp++; - Bp++; - } - } - return true; + return false; + } + Dp++; + Bp++; + } + } + return true; } template void Mat::simple_read (std::istream &st) { - double *Bp = (double*) data; - for (uint16 i=0;i> *Bp; - Bp++; - } - } + double *Bp = (double*) data; + for (uint16 i=0;i> *Bp; + Bp++; + } + } } @@ -401,113 +401,113 @@ class dMat; class dMat_interface { public: - dMat_interface ():ptr(NULL),M(0),N(0),row(0) - { } - double& operator[] (uint16 col) - { - assert(ptr); - assert (col < N); - assert (row < M); - return ptr[row*N + col]; - } - friend class dMat; + dMat_interface ():ptr(NULL),M(0),N(0),row(0) + { } + double& operator[] (uint16 col) + { + assert(ptr); + assert (col < N); + assert (row < M); + return ptr[row*N + col]; + } + friend class dMat; private: - uint16 M, N; - uint32 row; - double *ptr; + uint16 M, N; + uint32 row; + double *ptr; }; //dynamic matrix to pass arbitrary matrix to functions as arguments class dMat { public: - dMat(uint16 dim_m, uint16 dim_n) : dimM(0),dimN(0),data(NULL) - { - if (dim_m && dim_n) - resize(dim_m, dim_n); - } - dMat(uint16 dim_m, uint16 dim_n, double first, ...) : dimM(0),dimN(0),data(NULL) { + dMat(uint16 dim_m, uint16 dim_n) : dimM(0),dimN(0),data(NULL) + { + if (dim_m && dim_n) + resize(dim_m, dim_n); + } + dMat(uint16 dim_m, uint16 dim_n, double first, ...) : dimM(0),dimN(0),data(NULL) { va_list argp; - if (dim_m && dim_n) - { - resize(dim_m, dim_n); - va_start(argp, first); - data[0] = first; - for (uint16 i=1; i < dimM*dimN; i++) - data[i]=va_arg(argp, double); - va_end(argp); - } - } - dMat(const dMat &from) : dimM(0),dimN(0),data(NULL) - { - operator=(from); - } - dMat_interface& operator[] (uint16 row) { + if (dim_m && dim_n) + { + resize(dim_m, dim_n); + va_start(argp, first); + data[0] = first; + for (uint16 i=1; i < dimM*dimN; i++) + data[i]=va_arg(argp, double); + va_end(argp); + } + } + dMat(const dMat &from) : dimM(0),dimN(0),data(NULL) + { + operator=(from); + } + dMat_interface& operator[] (uint16 row) { //TODO: it seems unsafe for parallel read access to dMat.. - dmat_int.row = row; - return dmat_int; - } - void resize (uint16 dim_m, uint16 dim_n) - { - if (data) delete[] data; - data = new double[dim_m*dim_n]; - dmat_int.ptr = data; - dmat_int.M = dim_m; - dmat_int.N = dim_n; - dimM = dim_m; - dimN = dim_n; - zero(); - } - - void fill (double first, ...) { - va_list argp; - va_start(argp, first); - data[0] = first; - for (uint16 i=1; i < dimM*dimN; i++) - data[i]=va_arg(argp, double); - va_end(argp); - } - ~dMat() - { - if (data) delete[] data; - data = NULL; - } - - void zero () - { - memset(data,0,sizeof(double)*dimM*dimN); - } - template Mat toMat(); - template Vec toVec(uint16 col=0); - template void cpMat(Mat &mat); - template void cpVec(Vec &vec, uint16 col=0); - dMat& operator= (const dMat &from) - { - resize(from.dimM, from.dimN); - memcpy(this->data, from.data, sizeof(double)*dimM*dimN); - return *this; - } - - uint16 dM () - { - return dimM; - } - uint16 dN () - { - return dimN; - } + dmat_int.row = row; + return dmat_int; + } + void resize (uint16 dim_m, uint16 dim_n) + { + if (data) delete[] data; + data = new double[dim_m*dim_n]; + dmat_int.ptr = data; + dmat_int.M = dim_m; + dmat_int.N = dim_n; + dimM = dim_m; + dimN = dim_n; + zero(); + } + + void fill (double first, ...) { + va_list argp; + va_start(argp, first); + data[0] = first; + for (uint16 i=1; i < dimM*dimN; i++) + data[i]=va_arg(argp, double); + va_end(argp); + } + ~dMat() + { + if (data) delete[] data; + data = NULL; + } + + void zero () + { + memset(data,0,sizeof(double)*dimM*dimN); + } + template Mat toMat(); + template Vec toVec(uint16 col=0); + template void cpMat(Mat &mat); + template void cpVec(Vec &vec, uint16 col=0); + dMat& operator= (const dMat &from) + { + resize(from.dimM, from.dimN); + memcpy(this->data, from.data, sizeof(double)*dimM*dimN); + return *this; + } + + uint16 dM () + { + return dimM; + } + uint16 dN () + { + return dimN; + } double* ptr() { return data; } - friend std::ostream& operator<< (std::ostream& stream, dMat &obj); + friend std::ostream& operator<< (std::ostream& stream, dMat &obj); private: - uint16 dimM; - uint16 dimN; - double *data; - dMat_interface dmat_int; + uint16 dimM; + uint16 dimN; + double *data; + dMat_interface dmat_int; }; @@ -515,116 +515,116 @@ class dMat template Mat dMat::toMat() { - assert(M >= dimM && N >= dimN); //Mat >= dMat - Mat tmp; - for (uint16 i = 0; i < M; i++) - for (uint16 j = 0; j < N; j++) - tmp[i][j] = (*this)[i][j]; - return tmp; + assert(M >= dimM && N >= dimN); //Mat >= dMat + Mat tmp; + for (uint16 i = 0; i < M; i++) + for (uint16 j = 0; j < N; j++) + tmp[i][j] = (*this)[i][j]; + return tmp; } -template +template Vec dMat::toVec(uint16 col) { - assert(M >= dimM && dimN > col); //Vec >= dMat по M - Vec tmp; - for (uint16 i = 0; i < M; i++) - tmp[i] = (*this)[i][col]; - return tmp; + assert(M >= dimM && dimN > col); //Vec >= dMat по M + Vec tmp; + for (uint16 i = 0; i < M; i++) + tmp[i] = (*this)[i][col]; + return tmp; } -template +template void dMat::cpMat(Mat &mat) { - assert(dimM >= M && dimN >= N); //dMat >= Mat - for (uint16 i = 0; i < M; i++) - for (uint16 j = 0; j < N; j++) - (*this)[i][j] = mat[i][j]; + assert(dimM >= M && dimN >= N); //dMat >= Mat + for (uint16 i = 0; i < M; i++) + for (uint16 j = 0; j < N; j++) + (*this)[i][j] = mat[i][j]; } //функция копирует вектор в столбец col матрицы dMat // col - от нуля -template +template void dMat::cpVec(Vec &vec, uint16 col) { - assert(dimM >= M && dimN > col); - for (uint16 i = 0; i < M; i++) - (*this)[i][col] = vec[i]; + assert(dimM >= M && dimN > col); + for (uint16 i = 0; i < M; i++) + (*this)[i][col] = vec[i]; } //------------------------------------------------------------- -// MatSym +// MatSym //------------------------------------------------------------- template class MatSym { public: - MatSym () { - assert(dimM); - } - double* ptr() { - return data; - } - - const double* ptr() const { - return data; - } - - void zero() { - memset((void*) data, 0, sizeof(double)*getLength()); - } - + MatSym () { + assert(dimM); + } + double* ptr() { + return data; + } + + const double* ptr() const { + return data; + } + + void zero() { + memset((void*) data, 0, sizeof(double)*getLength()); + } + // indexing from 0 - double& comp (uint16 i, uint16 j) { - uint16 b; - if (i > j) { - b = j; - j = i; - i = b; - } - b = (2*dimM-(i-1))/2*i + (j-i); - b = dimM*i - (i-1)*i/2 + (j-i); - return data[b]; - } - - uint16 getLength () { - return dimM*(dimM+1)/2; - } + double& comp (uint16 i, uint16 j) { + uint16 b; + if (i > j) { + b = j; + j = i; + i = b; + } + b = (2*dimM-(i-1))/2*i + (j-i); + b = dimM*i - (i-1)*i/2 + (j-i); + return data[b]; + } + + uint16 getLength () { + return dimM*(dimM+1)/2; + } Mat toMat(); - void simple_read (std::istream &st); - bool compare (MatSym &B, double eps = 0.00001); - MatSym& operator+= (const MatSym &op); + void simple_read (std::istream &st); + bool compare (MatSym &B, double eps = 0.00001); + MatSym& operator+= (const MatSym &op); - double data[dimM*(dimM+1)/2]; + double data[dimM*(dimM+1)/2]; }; template void MatSym::simple_read (std::istream &st) { - double *Bp = (double*) data; - uint16 l = getLength(); - for (uint16 i=0;i> *Bp; - Bp++; - } + double *Bp = (double*) data; + uint16 l = getLength(); + for (uint16 i=0;i> *Bp; + Bp++; + } } template bool MatSym::compare (MatSym &B, double eps) { - double *Dp = (double*) data; - double *Bp = B.ptr(); - uint16 l = getLength(); - for (uint16 i=0;i eps) { - return false; - } - Dp++; - Bp++; - } - return true; + double *Dp = (double*) data; + double *Bp = B.ptr(); + uint16 l = getLength(); + for (uint16 i=0;i eps) { + return false; + } + Dp++; + Bp++; + } + return true; } template @@ -652,7 +652,7 @@ MatSym& MatSym::operator+= (const MatSym &op) { opp++; } - return *this; + return *this; } // [R] = coef * [B]^T*[D]*[B] @@ -661,59 +661,59 @@ MatSym& MatSym::operator+= (const MatSym &op) { // [D] - symmetric matrix (dimM x dimM) // [R] - symmetrix matrix (dimN x dimN) template -void matBTDBprod (Mat &B, MatSym &D, double coef, MatSym &R) +void matBTDBprod (Mat &B, MatSym &D, double coef, MatSym &R) { - Mat A; - A.zero(); - double *Ap = A.ptr(); - double *Bp = B.ptr(); - double *Dp = D.ptr(); - double *Rp = R.ptr(); - uint16 i,j,k; - - //A = BT*D - // BT*Du + BT*Dd - for (k=0;k A; + A.zero(); + double *Ap = A.ptr(); + double *Bp = B.ptr(); + double *Dp = D.ptr(); + double *Rp = R.ptr(); + uint16 i,j,k; + + //A = BT*D + // BT*Du + BT*Dd + for (k=0;k void matBTVprod(Mat &B, Vec &V, double coef, Vec &R) { #ifndef NLA3D_USE_BLAS - double *Bp = B.ptr(); - uint16 i,j; - for (i=0;i &B, Vec &V, double coef, Vec &R) template void matBVprod(Mat &B,Vec &V, double coef, Vec &R) { #ifndef NLA3D_USE_BLAS - double *Bp = B.ptr(); - uint16 i,j; - for (i=0;i &B,Vec &V, double coef, Vec &R) { // TODO: this this inefficient version // convert to regular matrix Mat mB = B.toMat(); - double *Bp = mB.ptr(); - uint16 i,j; - for (i=0;i void matABprod(Mat &A, Mat &B, const double coef, Mat &R) { #ifndef NLA3D_USE_BLAS - uint16 i,j,k; - double *Rp = R.ptr(); - double *Ap = A.ptr(); - double *Bp = B.ptr(); - for (i=0;i &A, Mat &B, const double coef, Mat< template void matATBprod(Mat &A, Mat &B, const double coef, Mat &R) { #ifndef NLA3D_USE_BLAS - uint16 i,j,k; - double *Rp = R.ptr(); - double *Ap = A.ptr(); - double *Bp = B.ptr(); - for (i=0;i 0 && _i <= nRows); - assert(_j > 0 && _j <= nColumns); + assert(_i > 0 && _i <= nRows); + assert(_j > 0 && _j <= nColumns); uint32 ind; - uint32 st = iofeir[_i-1] - 1; - uint32 en = iofeir[_i] - 1; + uint32 st = iofeir[_i-1] - 1; + uint32 en = iofeir[_i] - 1; if (st == en) return invalid; en--; - while(1) { - if (en - st == 1) { - if (columns[st] == _j) + while(1) { + if (en - st == 1) { + if (columns[st] == _j) return st; - if (columns[en] == _j) + if (columns[en] == _j) return en; - return invalid; - } - ind = (uint32) ((en+st)*0.5); - - if (columns[ind] == _j) + return invalid; + } + ind = (uint32) ((en+st)*0.5); + + if (columns[ind] == _j) return ind; - if (en == st) + if (en == st) return invalid; - if (columns[ind] > _j) { - en = ind; + if (columns[ind] > _j) { + en = ind; } else { - st = ind; + st = ind; } - } + } LOG(FATAL) << "What i'm doing here..?"; } @@ -159,25 +159,25 @@ void BaseSparseMatrix::printInternalData(std::ostream& out) { assert(values); assert(si->compressed); - out << "values = {"; - for (uint32 i = 0; i < si->numberOfValues; i++) { - out << values[i] << "\t"; + out << "values = {"; + for (uint32 i = 0; i < si->numberOfValues; i++) { + out << values[i] << "\t"; } - out << "}" << std::endl; + out << "}" << std::endl; - out << "columns = {"; - for (uint32 i = 0; i < si->nRows; i++) { - for (uint32 j = si->iofeir[i] - 1; j < si->iofeir[i + 1] - 1; j++) { - out << si->columns[j] << "\t"; + out << "columns = {"; + for (uint32 i = 0; i < si->nRows; i++) { + for (uint32 j = si->iofeir[i] - 1; j < si->iofeir[i + 1] - 1; j++) { + out << si->columns[j] << "\t"; } - } - out << "}" << std::endl; + } + out << "}" << std::endl; - out << "iofeir = {"; - for (uint32 i = 0; i <= si->nRows; i++) { - out << si->iofeir[i] << "\t"; + out << "iofeir = {"; + for (uint32 i = 0; i <= si->nRows; i++) { + out << si->iofeir[i] << "\t"; } - out << "}" << std::endl; + out << "}" << std::endl; } @@ -331,7 +331,7 @@ void BaseSparseMatrix::compress() { } else { si->compress(); } - + if (values) delete[] values; values = new double[si->numberOfValues]; zero(); @@ -344,7 +344,7 @@ void BaseSparseMatrix::setSparsityInfo(std::shared_ptr spar_info) // compressed = true state si = spar_info; // if we add SparsityInfo with alreade fixed number of entries we are ready to allocate values - if (si->compressed == true) + if (si->compressed == true) compress(); } @@ -381,14 +381,14 @@ void SparseMatrix::print(std::ostream& out) { assert(si->compressed); assert(values); - uint32 ind; - for (uint32 i = 1; i <= si->nRows; i++) { - out << "["; - for (uint32 j = 1; j <= si->nColumns; j++) { + uint32 ind; + for (uint32 i = 1; i <= si->nRows; i++) { + out << "["; + for (uint32 j = 1; j <= si->nColumns; j++) { out << value(i, j) << '\t'; } - out << "]" << std::endl; - } + out << "]" << std::endl; + } } @@ -396,7 +396,7 @@ double& SparseMatrix::operator()(uint32 _i, uint32 _j) { assert(values); assert(si); - uint32 index = si->getIndex(_i, _j); + uint32 index = si->getIndex(_i, _j); if (index == SparsityInfo::invalid) { LOG(FATAL) << "The position(" << _i << ", " << _j << ") is absent in the matrix"; } @@ -407,7 +407,7 @@ double SparseMatrix::value(uint32 _i, uint32 _j) const { assert(values); assert(si); - uint32 index = si->getIndex(_i, _j); + uint32 index = si->getIndex(_i, _j); if (index == SparsityInfo::invalid) { return 0.0; } @@ -453,15 +453,15 @@ void SparseSymMatrix::print(std::ostream& out) { assert(si->compressed); assert(values); - uint32 ind; - for (uint32 i = 1; i <= si->nRows; i++) { - // out << "["; - for (uint32 j = 1; j <= si->nColumns; j++) { + uint32 ind; + for (uint32 i = 1; i <= si->nRows; i++) { + // out << "["; + for (uint32 j = 1; j <= si->nColumns; j++) { out << value(i, j) << '\t'; } - // out << "]" << std::endl; + // out << "]" << std::endl; out << std::endl; - } + } } double& SparseSymMatrix::operator() (uint32 _i, uint32 _j) { @@ -469,9 +469,9 @@ double& SparseSymMatrix::operator() (uint32 _i, uint32 _j) { assert(si); // ensure that we work in upper triangle - if (_i > _j) std::swap(_i, _j); + if (_i > _j) std::swap(_i, _j); - uint32 index = si->getIndex(_i, _j); + uint32 index = si->getIndex(_i, _j); if (index == SparsityInfo::invalid) { LOG(FATAL) << "The position(" << _i << ", " << _j << ") is absent in the matrix"; } @@ -484,9 +484,9 @@ double SparseSymMatrix::value(uint32 _i, uint32 _j) const { assert(si); // ensure that we work in upper triangle - if (_i > _j) std::swap(_i, _j); + if (_i > _j) std::swap(_i, _j); - uint32 index = si->getIndex(_i, _j); + uint32 index = si->getIndex(_i, _j); if (index == SparsityInfo::invalid) { return 0.0; } @@ -502,7 +502,7 @@ void matBVprod(SparseSymMatrix &B, const dVec &V, const double coef, dVec &R) { assert(R.size() >= B.nRows()); // TODO: Try to use BLAS routines and measure speedup - const double eps = 1e-20; + const double eps = 1e-20; for (uint32 i = 1; i <= B.nRows(); i++) { // walk on lower triangle @@ -549,16 +549,16 @@ void matBTVprod(SparseMatrix &B, const dVec &V, const double coef, dVec &R) { assert(R.size() >= B.nColumns()); // TODO: Try to use BLAS routines and measure speedup - for (uint32 i = 1; i <= B.si->nRows; i++) { - if (B.si->iofeir[i] - B.si->iofeir[i-1] == 0) { - continue; + for (uint32 i = 1; i <= B.si->nRows; i++) { + if (B.si->iofeir[i] - B.si->iofeir[i-1] == 0) { + continue; } - uint32 st = B.si->iofeir[i-1] - 1; - uint32 en = B.si->iofeir[i] - 2; - for (uint32 j = st; j <= en; j++) + uint32 st = B.si->iofeir[i-1] - 1; + uint32 en = B.si->iofeir[i] - 2; + for (uint32 j = st; j <= en; j++) R[B.si->columns[j] - 1] += B.values[j] * V[i-1] * coef; - } + } } - + } // namespace math } // namespace nla3d diff --git a/src/lib/math/SparseMatrix.h b/src/lib/math/SparseMatrix.h index 86f7c78..171fe83 100755 --- a/src/lib/math/SparseMatrix.h +++ b/src/lib/math/SparseMatrix.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once @@ -89,7 +89,7 @@ class BaseSparseMatrix { BaseSparseMatrix(uint32 _nrows, uint32 _ncolumns, uint32 _max_in_row = 100); BaseSparseMatrix(std::shared_ptr spar_info); ~BaseSparseMatrix(); - + void reinit(uint32 _nrows, uint32 _ncols, const std::vector& entries); // perform compression procedure, after this new entries can't be added to matrix void compress(); @@ -136,10 +136,10 @@ class SparseMatrix : public BaseSparseMatrix { SparseMatrix(); SparseMatrix(uint32 nrows, uint32 ncolumns, uint32 max_in_row = 100); SparseMatrix(std::shared_ptr spar_info); - + void reinit(uint32 _nrows, uint32 _ncols, uint32 _max_in_row = 100); - // add non-zero entry to sparse matrix. This should be called before compress(). + // add non-zero entry to sparse matrix. This should be called before compress(). void addEntry(uint32 _i, uint32 _j); // add value to the _i, _j entry. This should be called after compress(). @@ -166,8 +166,8 @@ class SparseSymMatrix : public BaseSparseMatrix { SparseSymMatrix(std::shared_ptr spar_info); void reinit(uint32 _nrows, uint32 _max_in_row = 100); - - // add non-zero entry to sparse matrix. This should be called before compress(). + + // add non-zero entry to sparse matrix. This should be called before compress(). void addEntry(uint32 _i, uint32 _j); // add value to the _i, _j entry. This should be called after compress(). @@ -260,7 +260,7 @@ inline void SparseSymMatrix::addEntry(uint32 _i, uint32 _j) { assert(si); // ensure that we work in upper triangle - if (_i > _j) std::swap(_i, _j); + if (_i > _j) std::swap(_i, _j); si->addEntry(_i, _j); } @@ -268,6 +268,6 @@ inline void SparseSymMatrix::addValue(uint32 _i, uint32 _j, double value) { this->operator()(_i, _j) += value; } - + } // namespace math } // namespace nla3d diff --git a/src/lib/math/Vec.cpp b/src/lib/math/Vec.cpp index 69d9545..42d8f6c 100644 --- a/src/lib/math/Vec.cpp +++ b/src/lib/math/Vec.cpp @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #include "math/Vec.h" @@ -43,7 +43,7 @@ void dVec::reinit(dVec& _ref, uint32 _start, uint32 _size) { dVec::~dVec() { - clear(); + clear(); } uint32 dVec::size() const { @@ -101,7 +101,7 @@ double dVec::operator[](uint32 _n) const { dVec dVec::operator-() { assert(data); dVec p(size()); - for (uint32 i = 0; i < size(); i++) { + for (uint32 i = 0; i < size(); i++) { p[i] = -data[i]; } return p; @@ -115,7 +115,7 @@ dVec dVec::operator+(const dVec& op) { dVec p(size()); - for (uint32 i = 0; i < size(); i++) { + for (uint32 i = 0; i < size(); i++) { p[i] = data[i] + op.data[i]; } return p; @@ -130,7 +130,7 @@ dVec dVec::operator-(const dVec& op) { dVec p(size()); - for (uint32 i = 0; i < size(); i++) { + for (uint32 i = 0; i < size(); i++) { p[i] = data[i] - op.data[i]; } return p; @@ -142,7 +142,7 @@ dVec dVec::operator*(const double op) { dVec p(size()); - for (uint32 i = 0; i < size(); i++) { + for (uint32 i = 0; i < size(); i++) { p[i] = data[i] * op; } return p; @@ -153,7 +153,7 @@ dVec operator* (const double op1, const dVec& op2) { dVec p(op2.size()); - for (uint32 i = 0; i < op2.size(); i++) { + for (uint32 i = 0; i < op2.size(); i++) { p[i] = op2.data[i] * op1; } return p; @@ -165,7 +165,7 @@ dVec dVec::operator/(const double op) { dVec p(size()); - for (uint32 i = 0; i < size(); i++) { + for (uint32 i = 0; i < size(); i++) { p[i] = data[i] / op; } return p; @@ -178,7 +178,7 @@ dVec& dVec::operator+=(const dVec& op) { assert(op.data); assert(size() == op.size()); - for (uint32 i = 0; i < size(); i++) { + for (uint32 i = 0; i < size(); i++) { data[i] += op.data[i]; } return *this; @@ -190,7 +190,7 @@ dVec& dVec::operator-=(const dVec& op) { assert(op.data); assert(size() == op.size()); - for (uint32 i = 0; i < size(); i++) { + for (uint32 i = 0; i < size(); i++) { data[i] -= op.data[i]; } return *this; @@ -201,7 +201,7 @@ dVec& dVec::operator=(const dVec& op) { assert(op.data); assert(size() == op.size()); - for (uint32 i = 0; i < size(); i++) { + for (uint32 i = 0; i < size(); i++) { data[i] = op.data[i]; } return *this; diff --git a/src/lib/math/Vec.h b/src/lib/math/Vec.h index a91d1c8..e3d19e8 100755 --- a/src/lib/math/Vec.h +++ b/src/lib/math/Vec.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once #include "sys.h" @@ -12,184 +12,184 @@ namespace math { using namespace std; -template +template class Vec { public: - Vec() { - assert(dim); - memset(data,0,sizeof(double)*dim); //не очень красиво - } - Vec(double first, ...) { - assert(dim); - va_list argp; - va_start(argp, first); - data[0]=first; - for (uint16 i=1; i < dim; i++) data[i]=va_arg(argp, double); - va_end(argp); - } - ~Vec(void) { } - // TODO: this should be size_t instead of uint16 - double& operator[] (uint16 n); - const double operator[] (uint16 n) const; - void display (); - void zero() { - memset(data,0,sizeof(double)*dim); - } - double length (); - double qlength (); - Vec operator+ (const Vec &op); - Vec& operator+= (const Vec &op); - Vec operator- (); - Vec& operator= (const Vec &op); - Vec operator- (const Vec &op); - Vec operator* (const double op); - string toString (); - bool compare (Vec& V, double eps = 0.00005); - void simple_read(std::istream& st); - double operator* (const Vec &op); + Vec() { + assert(dim); + memset(data,0,sizeof(double)*dim); //не очень красиво + } + Vec(double first, ...) { + assert(dim); + va_list argp; + va_start(argp, first); + data[0]=first; + for (uint16 i=1; i < dim; i++) data[i]=va_arg(argp, double); + va_end(argp); + } + ~Vec(void) { } + // TODO: this should be size_t instead of uint16 + double& operator[] (uint16 n); + const double operator[] (uint16 n) const; + void display (); + void zero() { + memset(data,0,sizeof(double)*dim); + } + double length (); + double qlength (); + Vec operator+ (const Vec &op); + Vec& operator+= (const Vec &op); + Vec operator- (); + Vec& operator= (const Vec &op); + Vec operator- (const Vec &op); + Vec operator* (const double op); + string toString (); + bool compare (Vec& V, double eps = 0.00005); + void simple_read(std::istream& st); + double operator* (const Vec &op); template friend std::ostream& operator<< (std::ostream& stream,const Vec& obj); template friend Vec operator* (const double op1, const Vec &op2); - double* ptr (); + double* ptr (); private: - double data[dim]; + double data[dim]; }; //------operator[]------------------------------------------------------ template double& Vec::operator[] (uint16 n) { - assert(n < dim); - return data[n]; + assert(n < dim); + return data[n]; } //-------operator[]-const----------------------------------------------- template const double Vec::operator[] (uint16 n) const { - assert(n < dim); - return data[n]; + assert(n < dim); + return data[n]; } //--------display---------------------------------------------------- template void Vec::display () { - for (uint16 i = 0; i < dim; i++) { - std::cout << data[i]; - if (i < dim-1) std::cout << ", "; - } + for (uint16 i = 0; i < dim; i++) { + std::cout << data[i]; + if (i < dim-1) std::cout << ", "; + } } //-------operator+----------------------------------------------------- template Vec Vec::operator+(const Vec &op) { - Vec p; - for (uint16 i=0; i < dim; i++) p[i]=data[i]+op.data[i]; - return p; + Vec p; + for (uint16 i=0; i < dim; i++) p[i]=data[i]+op.data[i]; + return p; } //-------operator+=----------------------------------------------------- template Vec& Vec::operator+= (const Vec &op) { - for (uint16 i=0; i < dim; i++) this->data[i]+=op.data[i]; - return *this; + for (uint16 i=0; i < dim; i++) this->data[i]+=op.data[i]; + return *this; } //-------operator- ---------------------------------------------------- template Vec Vec::operator-() { - Vec p; - for (uint16 i=0; i < dim; i++) p[i]=-data[i]; - return p; + Vec p; + for (uint16 i=0; i < dim; i++) p[i]=-data[i]; + return p; } //------operator- ----------------------------------------------------- template Vec Vec::operator-(const Vec &op) { - Vec p; - for (uint16 i=0; i < dim; i++) p[i]=data[i]-op.data[i]; - return p; + Vec p; + for (uint16 i=0; i < dim; i++) p[i]=data[i]-op.data[i]; + return p; } //------operator*------------------------------------------------------ template Vec Vec::operator*(const double op) { - Vec p; - for (uint16 i=0; i < dim; i++) p[i]=data[i]*op; - return p; + Vec p; + for (uint16 i=0; i < dim; i++) p[i]=data[i]*op; + return p; } //------operator*------------------------------------------------------ template double Vec::operator*(const Vec &op) { - double p=0; - for (uint16 i=0; i < dim; i++) p+=data[i]*op.data[i]; - return p; + double p=0; + for (uint16 i=0; i < dim; i++) p+=data[i]*op.data[i]; + return p; } //---------qlength---------------------------------------------------- template double Vec::qlength() { - double p=0; - for (uint16 i=0; i < dim; i++) p+=data[i]*data[i]; - return p; + double p=0; + for (uint16 i=0; i < dim; i++) p+=data[i]*data[i]; + return p; } //---------length---------------------------------------------------- template double Vec::length() { - return sqrt(qlength()); + return sqrt(qlength()); } //---------operator*---------------------------------------------------- template Vec operator*(const double op1, const Vec &op2) { - Vec p; - for (uint16 i=0; i < dim1; i++) p[i]=op2.data[i]*op1; - return p; + Vec p; + for (uint16 i=0; i < dim1; i++) p[i]=op2.data[i]*op1; + return p; } //--------operator=------------------------------------------------------ template Vec& Vec::operator= (const Vec &op) { - memcpy(this->data, op.data, sizeof(double)*dim); - return *this; + memcpy(this->data, op.data, sizeof(double)*dim); + return *this; } //---------operator<<---------------------------------------------------- template std::ostream &operator<<(std::ostream &stream, const Vec &obj) { stream << obj.data[0]; - for (uint16 i = 1; i < dim1; i++) { - stream << " " << obj.data[i]; - } - return stream; + for (uint16 i = 1; i < dim1; i++) { + stream << " " << obj.data[i]; + } + return stream; } //------------------------------------------------------------- template string Vec::toString () { - string p; - char buff[100]; - for (uint16 i = 0; i < dim; i++) { - sprintf_s(buff,100,"%8.5e",this->data[i]); - p+=buff; - if (i < dim-1) p+= ", "; - } - return p; + string p; + char buff[100]; + for (uint16 i = 0; i < dim; i++) { + sprintf_s(buff,100,"%8.5e",this->data[i]); + p+=buff; + if (i < dim-1) p+= ", "; + } + return p; } //------------------------------------------------------- template double* Vec::ptr () { - return data; + return data; } //------------------------------------------------------- template bool Vec::compare (Vec& V, double eps) { - double *Dp = (double*) data; - double *Vp = V.ptr(); - for (uint16 j=0;j eps) { - return false; - } - Dp++; - Vp++; - } - return true; + double *Dp = (double*) data; + double *Vp = V.ptr(); + for (uint16 j=0;j eps) { + return false; + } + Dp++; + Vp++; + } + return true; } //------------------------------------------------------- template void Vec::simple_read(std::istream& st) { - double *Dp = data; - for (uint16 j=0;j> *Dp; - Dp++; - } + double *Dp = data; + for (uint16 j=0;j> *Dp; + Dp++; + } } diff --git a/src/lib/solidmech.cpp b/src/lib/solidmech.cpp index 51e8ba6..6063347 100644 --- a/src/lib/solidmech.cpp +++ b/src/lib/solidmech.cpp @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #include "solidmech.h" @@ -8,19 +8,19 @@ namespace nla3d { namespace solidmech { double J_C (const double* C) { - return sqrt( C[M_XX]*(C[M_YY]*C[M_ZZ]-C[M_YZ]*C[M_YZ]) + return sqrt( C[M_XX]*(C[M_YY]*C[M_ZZ]-C[M_YZ]*C[M_YZ]) - C[M_XY]*(C[M_XY]*C[M_ZZ]-C[M_YZ]*C[M_XZ]) + C[M_XZ]*(C[M_XY]*C[M_YZ]-C[M_YY]*C[M_XZ]) ); } void invC_C (const double* C, const double J, double* invC) { - double oo = 1.0/(J*J); - invC[M_XX] = oo*(C[M_YY]*C[M_ZZ]-C[M_YZ]*C[M_YZ]); - invC[M_XY] = oo*(C[M_XZ]*C[M_YZ]-C[M_XY]*C[M_ZZ]); - invC[M_XZ] = oo*(C[M_XY]*C[M_YZ]-C[M_XZ]*C[M_YY]); - invC[M_YY] = oo*(C[M_XX]*C[M_ZZ]-C[M_XZ]*C[M_XZ]); - invC[M_YZ] = oo*(C[M_XY]*C[M_XZ]-C[M_XX]*C[M_YZ]); - invC[M_ZZ] = oo*(C[M_XX]*C[M_YY]-C[M_XY]*C[M_XY]); + double oo = 1.0/(J*J); + invC[M_XX] = oo*(C[M_YY]*C[M_ZZ]-C[M_YZ]*C[M_YZ]); + invC[M_XY] = oo*(C[M_XZ]*C[M_YZ]-C[M_XY]*C[M_ZZ]); + invC[M_XZ] = oo*(C[M_XY]*C[M_YZ]-C[M_XZ]*C[M_YY]); + invC[M_YY] = oo*(C[M_XX]*C[M_ZZ]-C[M_XZ]*C[M_XZ]); + invC[M_YZ] = oo*(C[M_XY]*C[M_XZ]-C[M_XX]*C[M_YZ]); + invC[M_ZZ] = oo*(C[M_XX]*C[M_YY]-C[M_XY]*C[M_XY]); } @@ -35,10 +35,10 @@ void E_C (const double* C, double* E) { void IC_C (const double* C, double* IC) { - IC[0] = C[M_XX]+C[M_YY]+C[M_ZZ]; + IC[0] = C[M_XX]+C[M_YY]+C[M_ZZ]; //IC2 with star = 0.5 * IC[0]^2 - C:C - IC[1] = C[M_XX]*C[M_YY] + C[M_YY]*C[M_ZZ] + C[M_XX]*C[M_ZZ] - C[M_XY]*C[M_XY] - C[M_YZ]*C[M_YZ] - C[M_XZ]*C[M_XZ]; - IC[2] = J_C(C); + IC[1] = C[M_XX]*C[M_YY] + C[M_YY]*C[M_ZZ] + C[M_XX]*C[M_ZZ] - C[M_XY]*C[M_XY] - C[M_YZ]*C[M_YZ] - C[M_XZ]*C[M_XZ]; + IC[2] = J_C(C); } } // namespace nla3d diff --git a/src/lib/solidmech.h b/src/lib/solidmech.h index 7fe5307..9e1c8f9 100644 --- a/src/lib/solidmech.h +++ b/src/lib/solidmech.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once #include "sys.h" @@ -10,15 +10,15 @@ namespace solidmech { // labels for tensor components stored in 1-dim array enum tensorComponents { - M_XX = 0, - M_XY = 1, - M_XZ = 2, - M_YY = 3, - M_YZ = 4, - M_ZZ = 5 + M_XX = 0, + M_XY = 1, + M_XZ = 2, + M_YY = 3, + M_YZ = 4, + M_ZZ = 5 }; -// global order of tensor components in 1-dim array +// global order of tensor components in 1-dim array const tensorComponents defaultTensorComponents[] = {M_XX, M_XY, M_XZ, M_YY, M_YZ, M_ZZ}; const int LeviCivita[3][3][3] = { @@ -27,7 +27,7 @@ const int LeviCivita[3][3][3] = { {0,-1,0} }, { {0,0,-1}, {0,0,0}, - {1,0,0} }, + {1,0,0} }, { {0,1,0}, {-1,0,0}, {0,0,0} } }; diff --git a/src/lib/sys.cpp b/src/lib/sys.cpp index c9afa35..01b4385 100755 --- a/src/lib/sys.cpp +++ b/src/lib/sys.cpp @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #include "sys.h" @@ -12,68 +12,68 @@ namespace nla3d { static LogInitializer logInitializer; uint32 tick() { - return (uint32) clock(); + return (uint32) clock(); } int32 npow(int16 dig, uint16 power) { - int32 res=1; //TODO: if too big number? - for (uint16 i=0; i < power; i++) { - res *= dig; + int32 res=1; //TODO: if too big number? + for (uint16 i=0; i < power; i++) { + res *= dig; } - return res; + return res; } std::vector read_tokens(char *input) { - std::vector vec; - std::string tmp(""); - char delimeters[]="(),"; - char *p=input; - char *start=p; - while (*p) { - if (*p=='!') { + std::vector vec; + std::string tmp(""); + char delimeters[]="(),"; + char *p=input; + char *start=p; + while (*p) { + if (*p=='!') { break; } - bool isfind=false; - char* delp=delimeters; - while (*delp) { - if (*delp==*p) { - isfind=true; - break; - } - delp++; - } - if (isfind) { - tmp.assign(start, (int16) (p-start)); - del_spaces(tmp); - vec.push_back(tmp); - vec.push_back(std::string(delp,1)); - p++; - start=p; - continue; - } - p++; - } - tmp.assign(start, (int16) (p-start)); - del_spaces(tmp); - vec.push_back(tmp); - return vec; + bool isfind=false; + char* delp=delimeters; + while (*delp) { + if (*delp==*p) { + isfind=true; + break; + } + delp++; + } + if (isfind) { + tmp.assign(start, (int16) (p-start)); + del_spaces(tmp); + vec.push_back(tmp); + vec.push_back(std::string(delp,1)); + p++; + start=p; + continue; + } + p++; + } + tmp.assign(start, (int16) (p-start)); + del_spaces(tmp); + vec.push_back(tmp); + return vec; } void del_spaces (std::string &str) { - uint16 start = 0; - if (str.length()==0) return; - while (str[start]==' ' ) { - start++; - if (start == str.length()) { - str=" "; - return; - } - } - uint16 end = static_cast (str.length()-1); - while (str[end] ==' ') { + uint16 start = 0; + if (str.length()==0) return; + while (str[start]==' ' ) { + start++; + if (start == str.length()) { + str=" "; + return; + } + } + uint16 end = static_cast (str.length()-1); + while (str[end] ==' ') { end--; } - str=std::string(str,start, end-start+1); + str=std::string(str,start, end-start+1); } char* getCmdOption(char ** begin, char ** end, const std::string & option) { @@ -102,10 +102,10 @@ bool cmdOptionExists(char** begin, char** end, const std::string& option) { return std::find(begin, end, option) != end; } -//TODO: this functions only truncate file extension. +//TODO: this functions only truncate file extension. //But it was intended to leave only a file name (delete path and extension) std::string getFileNameFromPath(const std::string filename) { - std::string::const_reverse_iterator pivot = + std::string::const_reverse_iterator pivot = std::find( filename.rbegin(), filename.rend(), '.' ); return pivot == filename.rend() ? filename diff --git a/src/lib/sys.h b/src/lib/sys.h index cdabc52..5108d1e 100755 --- a/src/lib/sys.h +++ b/src/lib/sys.h @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #pragma once // before any include we need to define NOMINMAX to avoid redefining `max` with macros @@ -93,7 +93,7 @@ class LogInitializer { uint32 tick(); // template class to do convert from -// different types (mainly numerical) +// different types (mainly numerical) // to string. It should be just a wrapper // on C/C++ capabilities of conversation. // There are two options: @@ -144,7 +144,7 @@ uint16 str2dof (std::string dof_name); char* getCmdOption (char** begin, char** end, const std::string& option); bool cmdOptionExists (char** begin, char** end, const std::string& option); -std::vector getCmdManyOptions (char** begin, char** end, const std::string& option); +std::vector getCmdManyOptions (char** begin, char** end, const std::string& option); struct MatchPathSeparator { bool operator()( char ch ) const { diff --git a/src/main.cpp b/src/main.cpp index 1bd3758..e9b5d02 100755 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #include "sys.h" #include "FEStorage.h" @@ -13,7 +13,7 @@ using namespace nla3d; void prepeareProcessors (FEStorage& storage); -std::vector readRefCurveData (std::string fileRefCurve); +std::vector readRefCurveData (std::string fileRefCurve); double compareCurves (const std::vector& refCurve, const std::vector& curCurve); // Here is defaults values for command line options @@ -26,7 +26,7 @@ namespace options { bool useVtk = true; std::string modelFilename = ""; std::vector materialConstants; - std::string refCurveFilename = ""; + std::string refCurveFilename = ""; double curveCompareThreshold = 0.0; std::string reactionComponentName = ""; std::vector reactionDofs; @@ -95,7 +95,7 @@ bool parse_args (int argc, char* argv[]) { } else { options::reactionComponentName = vtmp[0]; for (size_t i = 1; i < vtmp.size(); i++) { - options::reactionDofs.push_back(Dof::label2dofType(vtmp[i])); + options::reactionDofs.push_back(Dof::label2dofType(vtmp[i])); } } } @@ -114,7 +114,7 @@ bool parse_args (int argc, char* argv[]) { options::rigidBodyDofs.push_back(Dof::UZ); } else { for (size_t i = 2; i < vtmp.size(); i++) { - options::rigidBodyDofs.push_back(Dof::label2dofType(vtmp[i])); + options::rigidBodyDofs.push_back(Dof::label2dofType(vtmp[i])); } } } @@ -284,8 +284,8 @@ int main (int argc, char* argv[]) { slaveDofLabels[i] = Dof::dofTypeLabels[i]; } - LOG(INFO) << "Rigid Body MPC collection: master node = " << mpc->masterNode << ", number of slaves = " << - mpc->slaveNodes.size() << ", slaves DoFs = " << slaveDofLabels; + LOG(INFO) << "Rigid Body MPC collection: master node = " << mpc->masterNode << ", number of slaves = " << + mpc->slaveNodes.size() << ", slaves DoFs = " << slaveDofLabels; } solver.solve(); @@ -307,7 +307,7 @@ int main (int argc, char* argv[]) { } LOG(INFO) << ss.str(); if (refCurve.size() > 0) { - curCurve = reactProc->getReactions(reactProc->dofs[0]); + curCurve = reactProc->getReactions(reactProc->dofs[0]); double _error = compareCurves (refCurve, curCurve); LOG(INFO) << "Error between reference loading curve and current is " << _error; if (options::curveCompareThreshold > 0.0) { @@ -332,7 +332,7 @@ double compareCurves (const std::vector& refCurve, const std::vector @@ -41,10 +41,10 @@ int main (int argc, char* argv[]) { } // Create an instance of FEStorage. - FEStorage storage; + FEStorage storage; // We have a deal with linear FE. Then it's ok to use linear solver (just one equilibrium iteration without // convergence controls) - LinearFESolver solver; + LinearFESolver solver; auto sind = storage.createNodes(md.nodesNumbers.size()); auto ind = md.nodesNumbers; @@ -96,8 +96,8 @@ int main (int argc, char* argv[]) { VtkProcessor* vtk = new VtkProcessor(&storage, "QUADTH"); solver.addPostProcessor(vtk); - solver.solve(); - + solver.solve(); + // Log all results about the model LOG(INFO) << "DoF solution:"; for (uint32 i = 1; i <= storage.nNodes(); i++) { @@ -115,7 +115,7 @@ int main (int argc, char* argv[]) { } } - return 0; + return 0; } diff --git a/test/QUADTH_transient.cpp b/test/QUADTH_transient.cpp index df8e152..21f2830 100644 --- a/test/QUADTH_transient.cpp +++ b/test/QUADTH_transient.cpp @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #include "sys.h" #include "FEStorage.h" @@ -13,29 +13,29 @@ using namespace nla3d; class ProbeProcessor : public PostProcessor { public: - ProbeProcessor(FEStorage *st); - virtual ~ProbeProcessor() { }; + ProbeProcessor(FEStorage *st); + virtual ~ProbeProcessor() { }; - virtual void pre (); - virtual void process(uint16 curLoadstep); - virtual void post(uint16 curLoadstep); + virtual void pre (); + virtual void process(uint16 curLoadstep); + virtual void post(uint16 curLoadstep); std::vector getMeasurments(); - - uint32 node = 0; - Dof::dofType dof = Dof::UNDEFINED; + + uint32 node = 0; + Dof::dofType dof = Dof::UNDEFINED; protected: - std::vector measurments; + std::vector measurments; }; ProbeProcessor::ProbeProcessor(FEStorage *st) : PostProcessor(st) { - name ="ProbeProcessor"; + name ="ProbeProcessor"; } void ProbeProcessor::pre() { - if (node == 0 || dof == Dof::UNDEFINED) { + if (node == 0 || dof == Dof::UNDEFINED) { LOG(FATAL) << "Can't work. No node number and/or Dof id. Processor name = " << name; return; } @@ -84,8 +84,8 @@ int main (int argc, char* argv[]) { } // Create an instance of FEStorage. - FEStorage storage; - LinearTransientFESolver solver; + FEStorage storage; + LinearTransientFESolver solver; auto sind = storage.createNodes(md.nodesNumbers.size()); auto ind = md.nodesNumbers; @@ -152,8 +152,8 @@ int main (int argc, char* argv[]) { probe->dof = Dof::TEMP; solver.addPostProcessor(probe); - solver.solve(); - + solver.solve(); + // Log all results about the model LOG(INFO) << " Probe Temperature:"; auto meas = probe->getMeasurments(); @@ -176,7 +176,7 @@ int main (int argc, char* argv[]) { CHECK(ave_fabs < 5.0e-2); } - return 0; + return 0; } diff --git a/test/TETRA0_test.cpp b/test/TETRA0_test.cpp index 0a4fdaf..3633dda 100644 --- a/test/TETRA0_test.cpp +++ b/test/TETRA0_test.cpp @@ -39,19 +39,19 @@ int main (int argc, char* argv[]) { MeshData md; if (!readCdbFile(cdb_filename, md)) { - LOG(FATAL) << "Can't read FE info from " << cdb_filename << "file. exiting.."; + LOG(FATAL) << "Can't read FE info from " << cdb_filename << "file. exiting.."; } md.compressNumbers(); - FEStorage storage; - LinearFESolver solver; + FEStorage storage; + LinearFESolver solver; - // add nodes - auto sind = storage.createNodes(md.nodesNumbers.size()); - auto ind = md.nodesNumbers; - for (uint32 i = 0; i < sind.size(); i++) { - storage.getNode(sind[i]).pos = md.nodesPos[i]; - } + // add nodes + auto sind = storage.createNodes(md.nodesNumbers.size()); + auto ind = md.nodesNumbers; + for (uint32 i = 0; i < sind.size(); i++) { + storage.getNode(sind[i]).pos = md.nodesPos[i]; + } ind = md.getCellsByAttribute("TYPE", 1); sind = storage.createElements(ind.size(), ElementType::TETRA0); @@ -79,7 +79,7 @@ int main (int argc, char* argv[]) { math::PARDISO_equationSolver eqSolver = math::PARDISO_equationSolver(); solver.attachEquationSolver(&eqSolver); #endif - solver.attachFEStorage (&storage); + solver.attachFEStorage (&storage); // VtkProcessor* vtk = new VtkProcessor (&storage, "tetra"); // solver.addPostProcessor(vtk); diff --git a/test/TimeControlTest.cpp b/test/TimeControlTest.cpp index 6d647b1..574c3dd 100644 --- a/test/TimeControlTest.cpp +++ b/test/TimeControlTest.cpp @@ -15,8 +15,8 @@ int main () { // Loadstep 5: dt = 15.0, currentTime = 75.0, numberOfEquilibrium = 5 // Loadstep 6: dt = 12.5, currentTime = 87.5, numberOfEquilibrium = 6 // Loadstep 7: dt = 12.5, currentTime = 100 , numberOfEquilibrium = 7 - - tc.setStartTime(0.0); + + tc.setStartTime(0.0); tc.setEndTime(100.0); CHECK(tc.getEndTime() == 100.0); CHECK(tc.getStartTime() == 0.0); diff --git a/test/eigen_mat_prod_test.cpp b/test/eigen_mat_prod_test.cpp index 85f4446..9f1b283 100644 --- a/test/eigen_mat_prod_test.cpp +++ b/test/eigen_mat_prod_test.cpp @@ -28,14 +28,14 @@ bool eigen_compare(const Eigen::MatrixBase& ref, const Eigen::MatrixBase } bool test_matBVprod () { - std::ifstream in; - const uint16 _N = 24; - const uint16 _M = 9; - char filename[100]; - Mat<_N,_M> B; - Vec<_M> V; - Vec<_N> Rf; - Vec<_N> R; + std::ifstream in; + const uint16 _N = 24; + const uint16 _M = 9; + char filename[100]; + Mat<_N,_M> B; + Vec<_M> V; + Vec<_N> Rf; + Vec<_N> R; Eigen::MatrixXd e_B; Eigen::VectorXd e_V; @@ -43,43 +43,43 @@ bool test_matBVprod () { sprintf_s(filename,100,"%s/matBVprod_%02d%02d",dir,_N,_M); in.open(filename); - for (uint16 gg = 1; gg <= tn; gg++) { - B.zero(); - V.zero(); - R.zero(); - Rf.zero(); + for (uint16 gg = 1; gg <= tn; gg++) { + B.zero(); + V.zero(); + R.zero(); + Rf.zero(); - B.simple_read(in); - V.simple_read(in); - Rf.simple_read(in); + B.simple_read(in); + V.simple_read(in); + Rf.simple_read(in); - matBVprod(B, V, 1.0, R); + matBVprod(B, V, 1.0, R); e_B = Eigen::Map > (B.ptr(), _N, _M); e_V = Eigen::Map (V.ptr(), _M, 1); e_R = e_B*e_V; - if (!R.compare(Rf, eps)) { - LOG(ERROR) << "test_matBVprod: n = " << gg; + if (!R.compare(Rf, eps)) { + LOG(ERROR) << "test_matBVprod: n = " << gg; exit(1); - } + } eigen_compare(Eigen::Map (Rf.ptr(),_N,1),e_R); - DLOG(DEBUG) << "test_matBVprod: case " << gg << " checked successfuly!"; - } + DLOG(DEBUG) << "test_matBVprod: case " << gg << " checked successfuly!"; + } in.close(); - return true; + return true; } bool test_matBTVprod () { - std::ifstream in; - const uint16 _N = 24; - const uint16 _M = 9; - char filename[100]; - Mat<_N,_M> B; - Vec<_N> V; - Vec<_M> Rf; - Vec<_M> R; + std::ifstream in; + const uint16 _N = 24; + const uint16 _M = 9; + char filename[100]; + Mat<_N,_M> B; + Vec<_N> V; + Vec<_M> Rf; + Vec<_M> R; Eigen::MatrixXd e_B; Eigen::VectorXd e_V; @@ -87,139 +87,139 @@ bool test_matBTVprod () { sprintf_s(filename,100, "%s/matBTVprod_%02d%02d", dir, _N, _M); in.open(filename); - for (uint16 gg = 1; gg <= tn; gg++) { - B.zero(); - V.zero(); - R.zero(); - Rf.zero(); + for (uint16 gg = 1; gg <= tn; gg++) { + B.zero(); + V.zero(); + R.zero(); + Rf.zero(); - B.simple_read(in); - V.simple_read(in); - Rf.simple_read(in); + B.simple_read(in); + V.simple_read(in); + Rf.simple_read(in); - matBTVprod(B, V, 1.0, R); + matBTVprod(B, V, 1.0, R); e_B = Eigen::Map > (B.ptr(), _N, _M); e_V = Eigen::Map (V.ptr(), _N, 1); e_R = e_B.transpose()*e_V; - if (!R.compare(Rf, eps)) { - LOG(ERROR) << "test_matBTVprod: n = " << gg; + if (!R.compare(Rf, eps)) { + LOG(ERROR) << "test_matBTVprod: n = " << gg; exit(1); - } + } eigen_compare(Eigen::Map (Rf.ptr(),_M,1), e_R); - DLOG(DEBUG) << "test_matBTVprod: case " << gg << " checked successfuly!"; - } + DLOG(DEBUG) << "test_matBTVprod: case " << gg << " checked successfuly!"; + } in.close(); - return true; + return true; } bool test_matABprod () { std::ifstream in; - const uint16 _N = 24; - const uint16 _M = 9; - const uint16 _M2= 12; - char filename[100]; - Mat<_N,_M> A; - Mat<_M,_M2> B; - Mat<_N,_M2> R; - Mat<_N,_M2> Rf; + const uint16 _N = 24; + const uint16 _M = 9; + const uint16 _M2= 12; + char filename[100]; + Mat<_N,_M> A; + Mat<_M,_M2> B; + Mat<_N,_M2> R; + Mat<_N,_M2> Rf; Eigen::MatrixXd e_A, e_B, e_R; sprintf_s(filename,100, "%s/matABprod_%02d%02d%02d", dir, _N, _M, _M2); in.open(filename); - for (uint16 gg = 1; gg <= tn; gg++) { - A.zero(); - B.zero(); - R.zero(); - Rf.zero(); + for (uint16 gg = 1; gg <= tn; gg++) { + A.zero(); + B.zero(); + R.zero(); + Rf.zero(); - A.simple_read(in); - B.simple_read(in); - Rf.simple_read(in); + A.simple_read(in); + B.simple_read(in); + Rf.simple_read(in); - matABprod(A, B, 1.0, R); + matABprod(A, B, 1.0, R); e_A = Eigen::Map > (A.ptr(), _N, _M); e_B = Eigen::Map > (B.ptr(), _M, _M2); e_R = e_A * e_B; - if (!R.compare(Rf, eps)) { - LOG(ERROR) << "test_matABprod: n = " << gg; + if (!R.compare(Rf, eps)) { + LOG(ERROR) << "test_matABprod: n = " << gg; exit(1); - } + } eigen_compare(Eigen::Map > (Rf.ptr(),_N, _M2), e_R); - DLOG(DEBUG) << "test_matABprod: case " << gg << " checked successfuly!"; - } + DLOG(DEBUG) << "test_matABprod: case " << gg << " checked successfuly!"; + } in.close(); - return true; + return true; } bool test_matATBprod () { std::ifstream in; - const uint16 _M = 24; - const uint16 _N = 9; - const uint16 _N2= 12; - char filename[100]; - Mat<_M,_N> A; - Mat<_M,_N2> B; - Mat<_N,_N2> R; - Mat<_N,_N2> Rf; + const uint16 _M = 24; + const uint16 _N = 9; + const uint16 _N2= 12; + char filename[100]; + Mat<_M,_N> A; + Mat<_M,_N2> B; + Mat<_N,_N2> R; + Mat<_N,_N2> Rf; Eigen::MatrixXd e_A, e_B, e_R; sprintf_s(filename,100, "%s/matATBprod_%02d%02d%02d", dir, _M, _N, _N2); in.open(filename); - for (uint16 gg = 1; gg <= tn; gg++) { - A.zero(); - B.zero(); - R.zero(); - Rf.zero(); + for (uint16 gg = 1; gg <= tn; gg++) { + A.zero(); + B.zero(); + R.zero(); + Rf.zero(); - A.simple_read(in); - B.simple_read(in); - Rf.simple_read(in); + A.simple_read(in); + B.simple_read(in); + Rf.simple_read(in); - matATBprod(A, B, 1.0, R); + matATBprod(A, B, 1.0, R); e_A = Eigen::Map > (A.ptr(), _M, _N); e_B = Eigen::Map > (B.ptr(), _M, _N2); e_R = e_A.transpose() * e_B; - if (!R.compare(Rf, eps)) { - LOG(ERROR) << "test_matATBprod: n = " << gg; + if (!R.compare(Rf, eps)) { + LOG(ERROR) << "test_matATBprod: n = " << gg; exit(1); - } + } eigen_compare(Eigen::Map > (Rf.ptr(),_N, _N2), e_R); - DLOG(DEBUG) << "test_matATBprod: case " << gg << " checked successfuly!"; - } + DLOG(DEBUG) << "test_matATBprod: case " << gg << " checked successfuly!"; + } in.close(); - return true; + return true; } bool test_matBTDBprod () { std::ifstream in; - const uint16 _M = 24; - const uint16 _N = 9; - char filename[100]; - Mat<_M,_N> B; - MatSym<_M> D; - MatSym<_N> R; - MatSym<_N> Rf; + const uint16 _M = 24; + const uint16 _N = 9; + char filename[100]; + Mat<_M,_N> B; + MatSym<_M> D; + MatSym<_N> R; + MatSym<_N> Rf; Eigen::MatrixXd e_D, e_B, e_R; sprintf_s(filename,100, "%s/matBTDBprod_%02d%02d", dir, _M, _N); in.open(filename); - for (uint16 gg = 1; gg <= tn; gg++) { - D.zero(); - B.zero(); - R.zero(); - Rf.zero(); + for (uint16 gg = 1; gg <= tn; gg++) { + D.zero(); + B.zero(); + R.zero(); + Rf.zero(); - B.simple_read(in); - D.simple_read(in); - Rf.simple_read(in); + B.simple_read(in); + D.simple_read(in); + Rf.simple_read(in); - matBTDBprod(B, D, 1.0, R); + matBTDBprod(B, D, 1.0, R); e_D.resize(_M,_M); @@ -247,14 +247,14 @@ bool test_matBTDBprod () { cout << "e_Rf.rows = " << e_Rf.rows() << ",e_Rf.cols = " << e_Rf.cols() << endl; cout << "e_R.rows = " << e_R.rows() << ",e_R.cols = " << e_R.cols() << endl; eigen_compare(e_Rf, e_R); - if (!R.compare(Rf, eps)) { - LOG(ERROR) << "test_matBTDBprod: n = " << gg; + if (!R.compare(Rf, eps)) { + LOG(ERROR) << "test_matBTDBprod: n = " << gg; exit(1); - } - DLOG(DEBUG) << "test_matBTDBprod: case " << gg << " checked successfuly!"; - } + } + DLOG(DEBUG) << "test_matBTDBprod: case " << gg << " checked successfuly!"; + } in.close(); - return true; + return true; } int main (int argc, char* argv[]) { @@ -269,9 +269,9 @@ int main (int argc, char* argv[]) { if (tmp) { tn = atoi(tmp); } - test_matBVprod(); - test_matBTVprod(); - test_matABprod(); - test_matATBprod(); - test_matBTDBprod(); + test_matBVprod(); + test_matBTVprod(); + test_matABprod(); + test_matATBprod(); + test_matBTDBprod(); } diff --git a/test/fereaders.cpp b/test/fereaders.cpp index 2a49649..a9a73b4 100644 --- a/test/fereaders.cpp +++ b/test/fereaders.cpp @@ -1,6 +1,6 @@ // 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 +// https://github.com/dmitryikh/nla3d #include "sys.h" #include "FEStorage.h" @@ -24,10 +24,10 @@ int main (int argc, char* argv[]) { // Create an instance of FEStorage. - FEStorage storage; + FEStorage storage; // We have a deal with linear FE. Then it's ok to use linear solver (just one equilibrium iteration without // convergence controls) - LinearFESolver solver; + LinearFESolver solver; MeshData mesh; if (!readCdbFile (cdb_filename, mesh)) { @@ -36,7 +36,7 @@ int main (int argc, char* argv[]) { LOG(INFO) << mesh.nodesNumbers.size() << " nodes were read"; for (size_t i = 0; i < mesh.nodesNumbers.size(); i++) { - LOG(INFO) << "Node " << mesh.nodesNumbers[i] << ": " << mesh.nodesPos[i][0] << " " + LOG(INFO) << "Node " << mesh.nodesNumbers[i] << ": " << mesh.nodesPos[i][0] << " " << mesh.nodesPos[i][1] << " " << mesh.nodesPos[i][2]; } @@ -47,8 +47,8 @@ int main (int argc, char* argv[]) { out << "Element " << mesh.cellNumbers[i] << ":"; for (size_t j = 0; j < mesh.cellNodes[i].size(); j++) { out << " " << mesh.cellNodes[i][j]; - } - LOG(INFO) << out.str(); + } + LOG(INFO) << out.str(); LOG(INFO) << " ETYPE:" << mesh.cellIntData["TYPE"][i] << " MAT:" << mesh.cellIntData["MAT"][i] << " REAL:" << mesh.cellIntData["REAL"][i] @@ -79,12 +79,12 @@ int main (int argc, char* argv[]) { // solver.attachEquationSolver(&eqSolver); // #endif // // FESolver should know FEStorage instance. Attach it. - // solver.attachFEStorage(&storage); + // solver.attachFEStorage(&storage); // VtkProcessor* vtk = new VtkProcessor(&storage, "QUADTH"); // solver.addPostProcessor(vtk); - // solver.solve(); - // + // solver.solve(); + // // // Log all results about the model // LOG(INFO) << "DoF solution:"; // for (uint32 i = 1; i <= storage.nNodes(); i++) { @@ -102,7 +102,7 @@ int main (int argc, char* argv[]) { // } // } - return 0; + return 0; } diff --git a/test/isoparametric_test.cpp b/test/isoparametric_test.cpp index 4955539..c7b40af 100644 --- a/test/isoparametric_test.cpp +++ b/test/isoparametric_test.cpp @@ -4,7 +4,7 @@ #include "elements/isoparametric.h" //TODO: not all results are autochecked -static const double th = 1.0e-8; +static const double th = 1.0e-8; using namespace std; using namespace nla3d; @@ -54,12 +54,12 @@ int main() { // arbitrary element nodes {-12.0, 1.0, 3.0}, // 5 {4.0, -5.0, -5.0}, // 6 -// arbitrary element nodes translated by +// arbitrary element nodes translated by {-12.0 + 9.0, +1.0 - 5.0, 3.0 - 3.0}, // 7 {+4.0 + 9.0, -5.0 -5.0, -5.0 - 3.0}}; // 8 // Create an instance of FEStorage. - FEStorage storage; + FEStorage storage; // Create and add nodes into FEStorage for (uint32 i = 1; i <= numberOfNodes; i++) { Node* no = new Node; @@ -98,7 +98,7 @@ int main() { // switch node numbers el = new LINEdummy; - el->getNodeNumber(0) = 6; + el->getNodeNumber(0) = 6; el->getNodeNumber(1) = 5; el->setIntegrationOrder(intOrder); storage.addElement(el); // 5 @@ -162,14 +162,14 @@ int main() { {4.0, -5.0, 0.0}, // 10 {7.0, 1.0, 0.0}, // 11 {-2.0, 3.0, 0.0}, // 12 -// arbitrary element nodes translated by +// arbitrary element nodes translated by {-12.0 + 9.0, 1.0 - 5.0, 0.0}, // 13 {4.0 + 9.0, -5.0 - 5.0, 0.0}, // 14 {7.0 + 9.0, 1.0 - 5.0, 0.0}, // 15 {-2.0 + 9.0, 3.0 - 5.0, 0.0}}; // 16 // Create an instance of FEStorage. - FEStorage storage; + FEStorage storage; // Create and add nodes into FEStorage for (uint32 i = 1; i <= numberOfNodes; i++) { Node* no = new Node; @@ -216,7 +216,7 @@ int main() { // switch node numbers el = new QUADdummy; - el->getNodeNumber(0) = 10; + el->getNodeNumber(0) = 10; el->getNodeNumber(1) = 11; el->getNodeNumber(2) = 12; el->getNodeNumber(3) = 9; @@ -228,7 +228,7 @@ int main() { cout << "Checking ElementIsoParamQUAD (-1,-1),(1,-1),(1,1),(-1,1) element.." << endl; el = dynamic_cast(&storage.getElement(1)); - for (uint32 np = 0; np < el->nOfIntPoints(); np++) { + for (uint32 np = 0; np < el->nOfIntPoints(); np++) { cout << "el->det[" << np << "] : " << el->det[np] << endl; CHECK_EQTH(el->det[np], 1.0, th); } @@ -238,7 +238,7 @@ int main() { cout << "Checking ElementIsoParamQUAD (-2,-1),(2,-1),(2,1),(-2,1) element.." << endl; el = dynamic_cast(&storage.getElement(2)); - for (uint32 np = 0; np < el->nOfIntPoints(); np++) { + for (uint32 np = 0; np < el->nOfIntPoints(); np++) { cout << "el->det[" << np << "] : " << el->det[np] << endl; CHECK_EQTH(el->det[np], 2.0, th); } @@ -252,13 +252,13 @@ int main() { QUADdummy* el3 = dynamic_cast(&storage.getElement(5)); assert(el1->getIntegrationOrder() == el2->getIntegrationOrder()); assert(el2->getIntegrationOrder() == el3->getIntegrationOrder()); - for (uint32 np = 0; np < el1->nOfIntPoints(); np++) { + for (uint32 np = 0; np < el1->nOfIntPoints(); np++) { cout << "el1->det[" << np << "] : " << el1->det[np] << endl; cout << "el2->det[" << np << "] : " << el2->det[np] << endl; cout << "el3->det[" << np << "] : " << el3->det[np] << endl; CHECK_EQTH(el1->det[np], el2->det[np], th); - for (uint16 no = 0; no < 4; no++) + for (uint16 no = 0; no < 4; no++) for (uint16 dim = 0; dim < 2; dim++) { // cout << "el1->NiXj[" << np << "][" << no << "][" << dim << "] " << el1->NiXj[np][no][dim] << endl; // cout << "el2->NiXj[" << np << "][" << no << "][" << dim << "] " << el2->NiXj[np][no][dim] << endl; @@ -272,7 +272,7 @@ int main() { CHECK_EQTH(el1->volume(), el2->volume(), th); CHECK_EQTH(el2->volume(), el3->volume(), th); } - + @@ -308,7 +308,7 @@ int main() { {4.0, -5.0, 4.0}, // 22 {7.0, 1.0, 4.0}, // 23 {-2.0, 3.0, 4.0}, // 24 -// arbitrary element nodes translated by +// arbitrary element nodes translated by {-12.0 + 9.0, 1.0 - 5.0, -7.0 + 2.0}, // 25 {4.0 + 9.0, -5.0 - 5.0, -7.0 + 2.0}, // 26 {7.0 + 9.0, 1.0 - 5.0, -7.0 + 2.0}, // 27 @@ -319,7 +319,7 @@ int main() { {-2.0 + 9.0, 3.0 - 5.0, 4.0 + 2.0}}; // 32 // Create an instance of FEStorage. - FEStorage storage; + FEStorage storage; // Create and add nodes into FEStorage for (uint32 i = 1; i <= numberOfNodes; i++) { Node* no = new Node; @@ -398,7 +398,7 @@ int main() { cout << "Checking ElementIsoParamHEXAHEDRON uni element.." << endl; el = dynamic_cast(&storage.getElement(1)); - for (uint32 np = 0; np < el->nOfIntPoints(); np++) { + for (uint32 np = 0; np < el->nOfIntPoints(); np++) { cout << "el->det[" << np << "] : " << el->det[np] << endl; CHECK_EQTH(el->det[np], 1.0, th); } @@ -408,7 +408,7 @@ int main() { cout << "Checking ElementIsoParamHEXAHEDRON stretched uni element.." << endl; el = dynamic_cast(&storage.getElement(2)); - for (uint32 np = 0; np < el->nOfIntPoints(); np++) { + for (uint32 np = 0; np < el->nOfIntPoints(); np++) { cout << "el->det[" << np << "] : " << el->det[np] << endl; CHECK_EQTH(el->det[np], 2.0, th); } @@ -422,13 +422,13 @@ int main() { HEXAHEDRONdummy* el3 = dynamic_cast(&storage.getElement(5)); assert(el1->getIntegrationOrder() == el2->getIntegrationOrder()); assert(el2->getIntegrationOrder() == el3->getIntegrationOrder()); - for (uint32 np = 0; np < el1->nOfIntPoints(); np++) { + for (uint32 np = 0; np < el1->nOfIntPoints(); np++) { cout << "el1->det[" << np << "] : " << el1->det[np] << endl; cout << "el2->det[" << np << "] : " << el2->det[np] << endl; cout << "el3->det[" << np << "] : " << el3->det[np] << endl; CHECK_EQTH(el1->det[np], el2->det[np], th); - for (uint16 no = 0; no < 8; no++) + for (uint16 no = 0; no < 8; no++) for (uint16 dim = 0; dim < 3; dim++) { // cout << "el1->NiXj[" << np << "][" << no << "][" << dim << "] " << el1->NiXj[np][no][dim] << endl; // cout << "el2->NiXj[" << np << "][" << no << "][" << dim << "] " << el2->NiXj[np][no][dim] << endl; diff --git a/test/mat_math.cpp b/test/mat_math.cpp index 6fa3331..0f9213c 100755 --- a/test/mat_math.cpp +++ b/test/mat_math.cpp @@ -9,7 +9,7 @@ using namespace nla3d::math; dMat test_func_dmat() { - return dMat((uint16) 3, (uint16) 3, (double) 1.1f, 4.4f, 7.7f, 5.5f, 2.2f, 8.8f, 6.6f, 9.9f, 3.3f); + return dMat((uint16) 3, (uint16) 3, (double) 1.1f, 4.4f, 7.7f, 5.5f, 2.2f, 8.8f, 6.6f, 9.9f, 3.3f); } bool compare_double(double d1, double d2) { @@ -17,80 +17,80 @@ bool compare_double(double d1, double d2) { return false; } return true; -} +} int main() { - Mat<3,3> mat(1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 9.0f); + Mat<3,3> mat(1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 9.0f); //test1 - Mat<3,3> mat2 = mat; - cout << "mat=mat2:" << endl; - cout << mat << endl; + Mat<3,3> mat2 = mat; + cout << "mat=mat2:" << endl; + cout << mat << endl; if (!mat.compare(mat2)) { LOG(ERROR) << "mat!=mat2"; exit(1); } //test2 - cout << "mat2.zero()" << endl; - mat2.zero(); - cout << mat2 << endl; + cout << "mat2.zero()" << endl; + mat2.zero(); + cout << mat2 << endl; if (!mat2.compare(Mat<3,3>(0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f))) { LOG(ERROR) << "mat2!=zeros"; exit(1); } - cout << "mat2.Identity()" << endl; - mat2.Identity(); - //mat=mat2; - cout << mat2 << endl; - cout << "mat*mat2" << endl; - cout << mat*mat2 << endl; - cout << "mat*mat" << endl; - cout << mat*mat << endl; - Mat<2,3> mat23(1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f); - cout << "mat23:" << endl << mat23; - cout << "mat23*mat" << endl << mat23*mat; - cout << "mat.transpose()" << endl << mat.transpose(); - Vec<3> vec(2.0f, 0.0f, 1.0f); - cout << "vec:" << endl << vec << endl; - cout << "mat*vec" << endl < cmat = mat; - cout << "cmat:" << endl << cmat << endl; - cout << "cmat[1][1]: " << cmat[1][1] << endl << endl; + cout << "mat2.Identity()" << endl; + mat2.Identity(); + //mat=mat2; + cout << mat2 << endl; + cout << "mat*mat2" << endl; + cout << mat*mat2 << endl; + cout << "mat*mat" << endl; + cout << mat*mat << endl; + Mat<2,3> mat23(1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f); + cout << "mat23:" << endl << mat23; + cout << "mat23*mat" << endl << mat23*mat; + cout << "mat.transpose()" << endl << mat.transpose(); + Vec<3> vec(2.0f, 0.0f, 1.0f); + cout << "vec:" << endl << vec << endl; + cout << "mat*vec" << endl < cmat = mat; + cout << "cmat:" << endl << cmat << endl; + cout << "cmat[1][1]: " << cmat[1][1] << endl << endl; //test Mat<3,3> - Mat<3,3> mat33(1.0, 2.0, 3.0, 43.0, 5.0, 6.0, 7.0, 8.0, 20.0); - cout << "mat33:" << endl << mat33 << endl; - cout << "det(mat33) = " << mat33.det() << " (right result = -657)" << endl; + Mat<3,3> mat33(1.0, 2.0, 3.0, 43.0, 5.0, 6.0, 7.0, 8.0, 20.0); + cout << "mat33:" << endl << mat33 << endl; + cout << "det(mat33) = " << mat33.det() << " (right result = -657)" << endl; if (!compare_double(mat33.det(), -657.0f)) { LOG(ERROR) << "mat33.det != -657.0"; exit(1); } - cout << "inv(mat33) = " << endl << mat33.inv(mat33.det()) << endl; + cout << "inv(mat33) = " << endl << mat33.inv(mat33.det()) << endl; Mat<3,3> inv_mat33(-0.07914764, 0.02435312, 0.00456621, 1.24505327, 0.00152207, -0.18721461, -0.47031963, -0.00913242, 0.12328767); cout << "answer inv_mat33 = " << endl << inv_mat33 << endl; if (!(mat33.inv(mat33.det())).compare(inv_mat33)) { LOG(ERROR) << "mat33.inv() is not correct"; exit(1); - } + } // test Mat<4,4> - Mat<4,4> mat44(1.0, 2.0, 3.0, 4.0, -5.0, 6.0, 7.0, 80.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 20.0); - cout << "mat44:" << endl << mat44 << endl; - cout << "det(mat44) = " << mat44.det() << " (right result = -320)" << endl; + Mat<4,4> mat44(1.0, 2.0, 3.0, 4.0, -5.0, 6.0, 7.0, 80.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 20.0); + cout << "mat44:" << endl << mat44 << endl; + cout << "det(mat44) = " << mat44.det() << " (right result = -320)" << endl; if (!compare_double(mat44.det(), -320.0f)) { LOG(ERROR) << "mat44.det != -320.0"; exit(1); } - cout << "inv(mat4) = " << endl << mat44.inv(mat44.det()) << endl; + cout << "inv(mat4) = " << endl << mat44.inv(mat44.det()) << endl; Mat<4,4> inv_mat44(9.50000000e-01, -1.00000000e-01, -2.65000000e+00, 1.80000000e+00, -3.15000000e+00, 2.00000000e-01, 5.30000000e+00, -3.35000000e+00, 1.95000000e+00, -1.00000000e-01, -2.15000000e+00, 1.30000000e+00, @@ -99,18 +99,18 @@ int main() if (!(mat44.inv(mat44.det())).compare(inv_mat44)) { LOG(ERROR) << "mat44.inv() is not correct"; exit(1); - } + } // test Mat<5,5> - Mat<5,5> mat55(1.0, 1.0, 3.0, 4.0, 5.0, 6.0, 7.0, 81.0, 9.0, 12.0, 11.0, 12.0, + Mat<5,5> mat55(1.0, 1.0, 3.0, 4.0, 5.0, 6.0, 7.0, 81.0, 9.0, 12.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 3.0, 18.0, 19.0, 20.0, 19.0, 22.0, 23.0, 24.0, 26.0); - cout << "mat55:" << endl << mat55 << endl; - cout << "det(mat55) = " << mat55.det() << " (right result = 1.001E4)" << endl; + cout << "mat55:" << endl << mat55 << endl; + cout << "det(mat55) = " << mat55.det() << " (right result = 1.001E4)" << endl; if (!compare_double(mat55.det(), 10005.0f)) { LOG(ERROR) << "mat55.det != 10005.0"; exit(1); } - cout << "inv(mat55) = " << endl << mat55.inv(mat55.det()) << endl; + cout << "inv(mat55) = " << endl << mat55.inv(mat55.det()) << endl; Mat<5,5> inv_mat55(-3.00649675e-01, -1.44927536e-02, -2.00239880e+00, 1.72613693e-01, 1.08695652e+00, -3.44827586e-02, -6.03330945e-18, 1.03448276e-01, -6.89655172e-02, @@ -123,48 +123,48 @@ int main() if (!(mat55.inv(mat55.det())).compare(inv_mat55)) { LOG(ERROR) << "mat55.inv() is not correct"; exit(1); - } - - Mat<3,6> MB(1.0f,2.0f,3.0f,4.0f,5.0f,6.0f,7.0f,8.0f,9.0f,10.0f,11.0f,12.0f,13.0f,14.0f,15.0f,16.0f,17.0f,18.0f); - Mat<3,3> D(11.0f,12.0f,13.0f,12.0f,22.0f,23.0f,13.0f,23.0f,33.0f); - cout << "MB:" << endl << MB; - cout << "D:" << endl << D; - cout << "MB.t*D*MB:" << endl << MB.transpose() * D* MB << endl; - Mat<6,3> tmp = MB.transpose()*D; - cout << "tmp = MB.t*D;tmp*MB" << endl << tmp*MB << endl; - Vec<3> ve(1.0f,2.0f,3.0f); - ve+=Vec<3>(10.0f, 20.0f, 30.0f); - cout << "ve:" << ve << endl; - - dMat dmat1(3,4); - cout << "dmat1:" << endl << dmat1 << endl; //выдает фигню. надо инициализировать массив нулями - dMat dmat2(3,3,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0); - cout << "dmat2:" << endl << dmat2 << endl; - cout << "mat33:" << endl << mat33 << endl;; - dmat1.cpMat<3,3>(mat33); - cout << "dmat1.cpMat<3,3>(mat33):" << endl << dmat1 << endl; - mat33 = dmat2.toMat<3,3>(); - cout << "mat33 = dmat2.toMat<3,3>()" << endl << mat33 << endl; - mat33 = test_func_dmat().toMat<3,3>(); - cout << "mat33 = test_func_dmat().toMat<3,3>()" << endl << mat33 << endl; - cout << "dmat1:" << endl << dmat1 << endl; - cout << "dmat2:" << endl << dmat2 << endl; - dmat1 = dmat2; - cout << "dmat1 = dmat2:" << endl << dmat1 << endl; + } + + Mat<3,6> MB(1.0f,2.0f,3.0f,4.0f,5.0f,6.0f,7.0f,8.0f,9.0f,10.0f,11.0f,12.0f,13.0f,14.0f,15.0f,16.0f,17.0f,18.0f); + Mat<3,3> D(11.0f,12.0f,13.0f,12.0f,22.0f,23.0f,13.0f,23.0f,33.0f); + cout << "MB:" << endl << MB; + cout << "D:" << endl << D; + cout << "MB.t*D*MB:" << endl << MB.transpose() * D* MB << endl; + Mat<6,3> tmp = MB.transpose()*D; + cout << "tmp = MB.t*D;tmp*MB" << endl << tmp*MB << endl; + Vec<3> ve(1.0f,2.0f,3.0f); + ve+=Vec<3>(10.0f, 20.0f, 30.0f); + cout << "ve:" << ve << endl; + + dMat dmat1(3,4); + cout << "dmat1:" << endl << dmat1 << endl; //выдает фигню. надо инициализировать массив нулями + dMat dmat2(3,3,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0); + cout << "dmat2:" << endl << dmat2 << endl; + cout << "mat33:" << endl << mat33 << endl;; + dmat1.cpMat<3,3>(mat33); + cout << "dmat1.cpMat<3,3>(mat33):" << endl << dmat1 << endl; + mat33 = dmat2.toMat<3,3>(); + cout << "mat33 = dmat2.toMat<3,3>()" << endl << mat33 << endl; + mat33 = test_func_dmat().toMat<3,3>(); + cout << "mat33 = test_func_dmat().toMat<3,3>()" << endl << mat33 << endl; + cout << "dmat1:" << endl << dmat1 << endl; + cout << "dmat2:" << endl << dmat2 << endl; + dmat1 = dmat2; + cout << "dmat1 = dmat2:" << endl << dmat1 << endl; // test some functions // Check matBVprod for symmetric matrices MatSym<3> sym3; sym3.zero(); - sym3.comp(0,0) = 5.0; - sym3.comp(0,1) = -35.0; - sym3.comp(0,2) = 12.0; + sym3.comp(0,0) = 5.0; + sym3.comp(0,1) = -35.0; + sym3.comp(0,2) = 12.0; - sym3.comp(1,1) = 15.0; - sym3.comp(1,2) = -3.0; - sym3.comp(2,2) = 14.0; + sym3.comp(1,1) = 15.0; + sym3.comp(1,2) = -3.0; + sym3.comp(2,2) = 14.0; - Vec<3> vec3(1.0f,2.0f,3.0f); + Vec<3> vec3(1.0f,2.0f,3.0f); Vec<3> res1; res1.zero(); @@ -181,13 +181,13 @@ int main() // check for operator+= for SymMat MatSym<3> sym3_2; sym3_2.zero(); - sym3_2.comp(0,0) = 7.0; - sym3_2.comp(0,1) = -25.0; - sym3_2.comp(0,2) = 19.0; + sym3_2.comp(0,0) = 7.0; + sym3_2.comp(0,1) = -25.0; + sym3_2.comp(0,2) = 19.0; - sym3_2.comp(1,1) = 0.0; - sym3_2.comp(1,2) = -7.0; - sym3_2.comp(2,2) = -3.3; + sym3_2.comp(1,1) = 0.0; + sym3_2.comp(1,2) = -7.0; + sym3_2.comp(2,2) = -3.3; sym3 += sym3_2; CHECK(sym3.comp(0,0) == 5.0 + 7.0); diff --git a/test/mat_prod_test.cpp b/test/mat_prod_test.cpp index ae59320..2a5ebe3 100755 --- a/test/mat_prod_test.cpp +++ b/test/mat_prod_test.cpp @@ -11,169 +11,169 @@ uint16 tn = 100; char* dir; bool test_matBVprod () { - std::ifstream in; - const uint16 _N = 24; - const uint16 _M = 9; - char filename[100]; - Mat<_N,_M> B; - Vec<_M> V; - Vec<_N> Rf; - Vec<_N> R; + std::ifstream in; + const uint16 _N = 24; + const uint16 _M = 9; + char filename[100]; + Mat<_N,_M> B; + Vec<_M> V; + Vec<_N> Rf; + Vec<_N> R; sprintf_s(filename,100,"%s/matBVprod_%02d%02d",dir,_N,_M); in.open(filename); - for (uint16 gg = 1; gg <= tn; gg++) { - B.zero(); - V.zero(); - R.zero(); - Rf.zero(); + for (uint16 gg = 1; gg <= tn; gg++) { + B.zero(); + V.zero(); + R.zero(); + Rf.zero(); - B.simple_read(in); - V.simple_read(in); - Rf.simple_read(in); + B.simple_read(in); + V.simple_read(in); + Rf.simple_read(in); - matBVprod(B, V, 1.0, R); + matBVprod(B, V, 1.0, R); - if (!R.compare(Rf, eps)) { - LOG(ERROR) << "test_matBVprod: n = " << gg; + if (!R.compare(Rf, eps)) { + LOG(ERROR) << "test_matBVprod: n = " << gg; exit(1); - } - DLOG(DEBUG) << "test_matBVprod: case " << gg << " checked successfuly!"; - } + } + DLOG(DEBUG) << "test_matBVprod: case " << gg << " checked successfuly!"; + } in.close(); - return true; + return true; } bool test_matBTVprod () { - std::ifstream in; - const uint16 _N = 24; - const uint16 _M = 9; - char filename[100]; - Mat<_N,_M> B; - Vec<_N> V; - Vec<_M> Rf; - Vec<_M> R; + std::ifstream in; + const uint16 _N = 24; + const uint16 _M = 9; + char filename[100]; + Mat<_N,_M> B; + Vec<_N> V; + Vec<_M> Rf; + Vec<_M> R; sprintf_s(filename,100, "%s/matBTVprod_%02d%02d", dir, _N, _M); in.open(filename); - for (uint16 gg = 1; gg <= tn; gg++) { - B.zero(); - V.zero(); - R.zero(); - Rf.zero(); + for (uint16 gg = 1; gg <= tn; gg++) { + B.zero(); + V.zero(); + R.zero(); + Rf.zero(); - B.simple_read(in); - V.simple_read(in); - Rf.simple_read(in); + B.simple_read(in); + V.simple_read(in); + Rf.simple_read(in); - matBTVprod(B, V, 1.0, R); + matBTVprod(B, V, 1.0, R); - if (!R.compare(Rf, eps)) { - LOG(ERROR) << "test_matBTVprod: n = " << gg; + if (!R.compare(Rf, eps)) { + LOG(ERROR) << "test_matBTVprod: n = " << gg; exit(1); - } - DLOG(DEBUG) << "test_matBTVprod: case " << gg << " checked successfuly!"; - } + } + DLOG(DEBUG) << "test_matBTVprod: case " << gg << " checked successfuly!"; + } in.close(); - return true; + return true; } bool test_matABprod () { std::ifstream in; - const uint16 _N = 24; - const uint16 _M = 9; - const uint16 _M2= 12; - char filename[100]; - Mat<_N,_M> A; - Mat<_M,_M2> B; - Mat<_N,_M2> R; - Mat<_N,_M2> Rf; + const uint16 _N = 24; + const uint16 _M = 9; + const uint16 _M2= 12; + char filename[100]; + Mat<_N,_M> A; + Mat<_M,_M2> B; + Mat<_N,_M2> R; + Mat<_N,_M2> Rf; sprintf_s(filename,100, "%s/matABprod_%02d%02d%02d", dir, _N, _M, _M2); in.open(filename); - for (uint16 gg = 1; gg <= tn; gg++) { - A.zero(); - B.zero(); - R.zero(); - Rf.zero(); + for (uint16 gg = 1; gg <= tn; gg++) { + A.zero(); + B.zero(); + R.zero(); + Rf.zero(); - A.simple_read(in); - B.simple_read(in); - Rf.simple_read(in); + A.simple_read(in); + B.simple_read(in); + Rf.simple_read(in); - matABprod(A, B, 1.0, R); + matABprod(A, B, 1.0, R); - if (!R.compare(Rf, eps)) { - LOG(ERROR) << "test_matABprod: n = " << gg; + if (!R.compare(Rf, eps)) { + LOG(ERROR) << "test_matABprod: n = " << gg; exit(1); - } - DLOG(DEBUG) << "test_matABprod: case " << gg << " checked successfuly!"; - } + } + DLOG(DEBUG) << "test_matABprod: case " << gg << " checked successfuly!"; + } in.close(); - return true; + return true; } bool test_matATBprod () { std::ifstream in; - const uint16 _M = 24; - const uint16 _N = 9; - const uint16 _N2= 12; - char filename[100]; - Mat<_M,_N> A; - Mat<_M,_N2> B; - Mat<_N,_N2> R; - Mat<_N,_N2> Rf; + const uint16 _M = 24; + const uint16 _N = 9; + const uint16 _N2= 12; + char filename[100]; + Mat<_M,_N> A; + Mat<_M,_N2> B; + Mat<_N,_N2> R; + Mat<_N,_N2> Rf; sprintf_s(filename,100, "%s/matATBprod_%02d%02d%02d", dir, _M, _N, _N2); in.open(filename); - for (uint16 gg = 1; gg <= tn; gg++) { - A.zero(); - B.zero(); - R.zero(); - Rf.zero(); - - A.simple_read(in); - B.simple_read(in); - Rf.simple_read(in); - - matATBprod(A, B, 1.0, R); - if (!R.compare(Rf, eps)) { - LOG(ERROR) << "test_matATBprod: n = " << gg; + for (uint16 gg = 1; gg <= tn; gg++) { + A.zero(); + B.zero(); + R.zero(); + Rf.zero(); + + A.simple_read(in); + B.simple_read(in); + Rf.simple_read(in); + + matATBprod(A, B, 1.0, R); + if (!R.compare(Rf, eps)) { + LOG(ERROR) << "test_matATBprod: n = " << gg; exit(1); - } - DLOG(DEBUG) << "test_matATBprod: case " << gg << " checked successfuly!"; - } + } + DLOG(DEBUG) << "test_matATBprod: case " << gg << " checked successfuly!"; + } in.close(); - return true; + return true; } bool test_matBTDBprod () { std::ifstream in; - const uint16 _M = 24; - const uint16 _N = 9; - char filename[100]; - Mat<_M,_N> B; - MatSym<_M> D; - MatSym<_N> R; - MatSym<_N> Rf; + const uint16 _M = 24; + const uint16 _N = 9; + char filename[100]; + Mat<_M,_N> B; + MatSym<_M> D; + MatSym<_N> R; + MatSym<_N> Rf; sprintf_s(filename,100, "%s/matBTDBprod_%02d%02d", dir, _M, _N); in.open(filename); - for (uint16 gg = 1; gg <= tn; gg++) { - D.zero(); - B.zero(); - R.zero(); - Rf.zero(); + for (uint16 gg = 1; gg <= tn; gg++) { + D.zero(); + B.zero(); + R.zero(); + Rf.zero(); - B.simple_read(in); - D.simple_read(in); - Rf.simple_read(in); + B.simple_read(in); + D.simple_read(in); + Rf.simple_read(in); - matBTDBprod(B, D, 1.0, R); + matBTDBprod(B, D, 1.0, R); - if (!R.compare(Rf, eps)) { - LOG(ERROR) << "test_matBTDBprod: n = " << gg; + if (!R.compare(Rf, eps)) { + LOG(ERROR) << "test_matBTDBprod: n = " << gg; exit(1); - } - DLOG(DEBUG) << "test_matBTDBprod: case " << gg << " checked successfuly!"; - } + } + DLOG(DEBUG) << "test_matBTDBprod: case " << gg << " checked successfuly!"; + } in.close(); - return true; + return true; } int main (int argc, char* argv[]) { @@ -188,9 +188,9 @@ int main (int argc, char* argv[]) { if (tmp) { tn = atoi(tmp); } - test_matBVprod(); - test_matBTVprod(); - test_matABprod(); - test_matATBprod(); - test_matBTDBprod(); + test_matBVprod(); + test_matBTVprod(); + test_matABprod(); + test_matATBprod(); + test_matBTDBprod(); } diff --git a/test/mat_sparse_pardiso.cpp b/test/mat_sparse_pardiso.cpp index 2be9d1a..43b0e80 100755 --- a/test/mat_sparse_pardiso.cpp +++ b/test/mat_sparse_pardiso.cpp @@ -7,154 +7,154 @@ using namespace nla3d::math; //TODO: undone! (see commented lines below) int main() { - int n = 8; - SparseSymmetricMatrix spmat(8); - spmat.startTraining(); - //// 7.0 0.0 1.0 0.0 0.0 2.0 7.0 0.0 - //// 0.0 -4.0 8.0 0.0 2.0 0.0 0.0 0.0 - //// 1.0 8.0 1.0 0.0 0.0 0.0 0.0 5.0 - //// 0.0 0.0 0.0 7.0 0.0 0.0 9.0 0.0 - //// 0.0 2.0 0.0 0.0 5.0 1.0 5.0 0.0 - //// 2.0 0.0 0.0 0.0 1.0 -1.0 0.0 5.0 - //// 7.0 0.0 0.0 9.0 5.0 0.0 11.0 0.0 - //// 0.0 0.0 5.0 0.0 0.0 5.0 0.0 5.0 - spmat.addValue(1,1,7.0); - spmat.addValue(1,3,1.0); - spmat.addValue(1,6,2.0); - spmat.addValue(1,7,7.0); - spmat.addValue(2,5,2.0); - spmat.addValue(2,3,8.0); - spmat.addValue(2,2,-4.0); - spmat.addValue(3,3,1.0); - spmat.addValue(3,8,5.0); - spmat.addValue(4,7,9.0); - spmat.addValue(4,4,7.0); - spmat.addValue(5,5,5.0); - spmat.addValue(5,6,1.0); - spmat.addValue(5,7,5.0); - spmat.addValue(6,8,5.0); - spmat.addValue(6,6,-1.0); - spmat.addValue(7,7,11.0); - spmat.addValue(8,8,5.0); - spmat.stopTraining(); - //spmat.print_data(); // works correctly - double b[] = {1.0, 2.0, 4.0, 6.0, 5.0, 3.0, 7.0, 3.0}; - cout << "SUR*vec: " << endl; - for (int32 i=0; i < 8; i++) - cout << spmat.mult_vec_i(b,i+1) << endl; + int n = 8; + SparseSymmetricMatrix spmat(8); + spmat.startTraining(); + //// 7.0 0.0 1.0 0.0 0.0 2.0 7.0 0.0 + //// 0.0 -4.0 8.0 0.0 2.0 0.0 0.0 0.0 + //// 1.0 8.0 1.0 0.0 0.0 0.0 0.0 5.0 + //// 0.0 0.0 0.0 7.0 0.0 0.0 9.0 0.0 + //// 0.0 2.0 0.0 0.0 5.0 1.0 5.0 0.0 + //// 2.0 0.0 0.0 0.0 1.0 -1.0 0.0 5.0 + //// 7.0 0.0 0.0 9.0 5.0 0.0 11.0 0.0 + //// 0.0 0.0 5.0 0.0 0.0 5.0 0.0 5.0 + spmat.addValue(1,1,7.0); + spmat.addValue(1,3,1.0); + spmat.addValue(1,6,2.0); + spmat.addValue(1,7,7.0); + spmat.addValue(2,5,2.0); + spmat.addValue(2,3,8.0); + spmat.addValue(2,2,-4.0); + spmat.addValue(3,3,1.0); + spmat.addValue(3,8,5.0); + spmat.addValue(4,7,9.0); + spmat.addValue(4,4,7.0); + spmat.addValue(5,5,5.0); + spmat.addValue(5,6,1.0); + spmat.addValue(5,7,5.0); + spmat.addValue(6,8,5.0); + spmat.addValue(6,6,-1.0); + spmat.addValue(7,7,11.0); + spmat.addValue(8,8,5.0); + spmat.stopTraining(); + //spmat.print_data(); // works correctly + double b[] = {1.0, 2.0, 4.0, 6.0, 5.0, 3.0, 7.0, 3.0}; + cout << "SUR*vec: " << endl; + for (int32 i=0; i < 8; i++) + cout << spmat.mult_vec_i(b,i+1) << endl; - SparseMatrix spmat2(8,8); - spmat2.startTraining(); - //// 7.0 0.0 1.0 0.0 0.0 2.0 7.0 0.0 - //// 0.0 -4.0 8.0 0.0 2.0 0.0 0.0 0.0 - //// 1.0 8.0 1.0 0.0 0.0 0.0 0.0 5.0 - //// 0.0 0.0 0.0 7.0 0.0 0.0 9.0 0.0 - //// 0.0 2.0 0.0 0.0 5.0 1.0 5.0 0.0 - //// 2.0 0.0 0.0 0.0 1.0 -1.0 0.0 5.0 - //// 7.0 0.0 0.0 9.0 5.0 0.0 11.0 0.0 - //// 0.0 0.0 5.0 0.0 0.0 5.0 0.0 5.0 - spmat2.addValue(1,1,7.0); - spmat2.addValue(3,1,1.0); - spmat2.addValue(1,6,2.0); - spmat2.addValue(1,7,7.0); - spmat2.addValue(2,5,2.0); - spmat2.addValue(2,3,8.0); - spmat2.addValue(2,2,-4.0); - spmat2.addValue(3,3,1.0); - spmat2.addValue(3,8,5.0); - spmat2.addValue(4,7,9.0); - spmat2.addValue(4,4,7.0); - spmat2.addValue(5,5,5.0); - spmat2.addValue(5,6,1.0); - spmat2.addValue(5,7,5.0); - spmat2.addValue(6,8,5.0); - spmat2.addValue(6,6,-1.0); - spmat2.addValue(7,7,11.0); - spmat2.addValue(8,8,5.0); - spmat2.addValue(7,1,7.0); - spmat2.stopTraining(); - cout << "R*vec: " << endl; - for (int32 i=0; i < 8; i++) - cout << spmat2.mult_vec_i(b,i+1) << endl; - cout << "R^T*vec: " << endl; - for (int32 i=0; i < 8; i++) - cout << spmat2.transpose_mult_vec_i(b,i+1) << endl; - //double x[8]; - //int nrhs = 1; - //int mtype = -2; //real symmetric matrix - //void *pt[64]; /*Internal solver memory pointer pt,*/ - //int iparm[64]; - //int maxfct, mnum, phase, error, msglvl; - //int i; - //for (i=0;i<64;i++) iparm[i]=0; - //iparm[0] = 1; //no solver default - //iparm[1] = 2; //fill-in reordering from meris - //iparm[2] = MKL_Get_Max_Threads(); - //iparm[3] = 0; //no iterative-direct algorithm - //iparm[4] = 0; //no user fill-in reducing permutation - //iparm[5] = 0; //write solution into x - //iparm[6] = 16; //default logical fortran unit number for output - //iparm[7] = 2; //max numbers of iterative refinement steps - //iparm[9] = 13; //pertrub the pivor elements with 1e-13 - //iparm[10] = 1; //use nonsymmetric permutation and scaling MPS - //iparm[13]=0; //output: number of perturbed pivots - //iparm[17]=-1; //output: number of nonzeros in the factor LU - //iparm[18]=-1; //output: MFLOPS for LU factorization - //iparm[19] = 0; //output: number of CG Iterations + SparseMatrix spmat2(8,8); + spmat2.startTraining(); + //// 7.0 0.0 1.0 0.0 0.0 2.0 7.0 0.0 + //// 0.0 -4.0 8.0 0.0 2.0 0.0 0.0 0.0 + //// 1.0 8.0 1.0 0.0 0.0 0.0 0.0 5.0 + //// 0.0 0.0 0.0 7.0 0.0 0.0 9.0 0.0 + //// 0.0 2.0 0.0 0.0 5.0 1.0 5.0 0.0 + //// 2.0 0.0 0.0 0.0 1.0 -1.0 0.0 5.0 + //// 7.0 0.0 0.0 9.0 5.0 0.0 11.0 0.0 + //// 0.0 0.0 5.0 0.0 0.0 5.0 0.0 5.0 + spmat2.addValue(1,1,7.0); + spmat2.addValue(3,1,1.0); + spmat2.addValue(1,6,2.0); + spmat2.addValue(1,7,7.0); + spmat2.addValue(2,5,2.0); + spmat2.addValue(2,3,8.0); + spmat2.addValue(2,2,-4.0); + spmat2.addValue(3,3,1.0); + spmat2.addValue(3,8,5.0); + spmat2.addValue(4,7,9.0); + spmat2.addValue(4,4,7.0); + spmat2.addValue(5,5,5.0); + spmat2.addValue(5,6,1.0); + spmat2.addValue(5,7,5.0); + spmat2.addValue(6,8,5.0); + spmat2.addValue(6,6,-1.0); + spmat2.addValue(7,7,11.0); + spmat2.addValue(8,8,5.0); + spmat2.addValue(7,1,7.0); + spmat2.stopTraining(); + cout << "R*vec: " << endl; + for (int32 i=0; i < 8; i++) + cout << spmat2.mult_vec_i(b,i+1) << endl; + cout << "R^T*vec: " << endl; + for (int32 i=0; i < 8; i++) + cout << spmat2.transpose_mult_vec_i(b,i+1) << endl; + //double x[8]; + //int nrhs = 1; + //int mtype = -2; //real symmetric matrix + //void *pt[64]; /*Internal solver memory pointer pt,*/ + //int iparm[64]; + //int maxfct, mnum, phase, error, msglvl; + //int i; + //for (i=0;i<64;i++) iparm[i]=0; + //iparm[0] = 1; //no solver default + //iparm[1] = 2; //fill-in reordering from meris + //iparm[2] = MKL_Get_Max_Threads(); + //iparm[3] = 0; //no iterative-direct algorithm + //iparm[4] = 0; //no user fill-in reducing permutation + //iparm[5] = 0; //write solution into x + //iparm[6] = 16; //default logical fortran unit number for output + //iparm[7] = 2; //max numbers of iterative refinement steps + //iparm[9] = 13; //pertrub the pivor elements with 1e-13 + //iparm[10] = 1; //use nonsymmetric permutation and scaling MPS + //iparm[13]=0; //output: number of perturbed pivots + //iparm[17]=-1; //output: number of nonzeros in the factor LU + //iparm[18]=-1; //output: MFLOPS for LU factorization + //iparm[19] = 0; //output: number of CG Iterations - //maxfct = 1; //maximum number of numerical factorizations - //mnum = 1; //wich factorization to use - //msglvl = 0; //dont print statistical information in file - //error = 0; //init error flag - //for (i=0;i<64;i++) pt[i]=0; + //maxfct = 1; //maximum number of numerical factorizations + //mnum = 1; //wich factorization to use + //msglvl = 0; //dont print statistical information in file + //error = 0; //init error flag + //for (i=0;i<64;i++) pt[i]=0; - ////Reordering and Symbolic Factorization. This step also allocate - //// all memory that is necessary for the factorization - //phase = 11; - //PARDISO(pt, &maxfct, &mnum, &mtype,&phase, - // &n, spmat.get_values_array(), (int*) spmat.get_iofeir_array(), (int*) spmat.get_columns_array(), NULL, &nrhs, iparm, &msglvl, NULL, NULL, &error); - //if (error != 0) - //{ - // printf("\nERROR during symbolic factorization: %d", error); - // exit(1); - //} + ////Reordering and Symbolic Factorization. This step also allocate + //// all memory that is necessary for the factorization + //phase = 11; + //PARDISO(pt, &maxfct, &mnum, &mtype,&phase, + // &n, spmat.get_values_array(), (int*) spmat.get_iofeir_array(), (int*) spmat.get_columns_array(), NULL, &nrhs, iparm, &msglvl, NULL, NULL, &error); + //if (error != 0) + //{ + // printf("\nERROR during symbolic factorization: %d", error); + // exit(1); + //} - //printf("\nReordering completed..."); - //printf("\nNumber of nonzeros in factors=%d", iparm[17]); - //printf("\nNumber of factorization MFLOPS=%d", iparm[18]); + //printf("\nReordering completed..."); + //printf("\nNumber of nonzeros in factors=%d", iparm[17]); + //printf("\nNumber of factorization MFLOPS=%d", iparm[18]); - //// Numerical factorization + //// Numerical factorization - //phase = 22; - //PARDISO(pt, &maxfct, &mnum, &mtype,&phase, - // &n, spmat.get_values_array(), (int*) spmat.get_iofeir_array(), (int*) spmat.get_columns_array(), NULL, &nrhs, iparm, &msglvl, NULL, NULL, &error); - //if (error != 0) - //{ - // printf("\nERROR during numerical factorization: %d", error); - // exit(2); - //} - //printf("\nFactorization completed..."); + //phase = 22; + //PARDISO(pt, &maxfct, &mnum, &mtype,&phase, + // &n, spmat.get_values_array(), (int*) spmat.get_iofeir_array(), (int*) spmat.get_columns_array(), NULL, &nrhs, iparm, &msglvl, NULL, NULL, &error); + //if (error != 0) + //{ + // printf("\nERROR during numerical factorization: %d", error); + // exit(2); + //} + //printf("\nFactorization completed..."); - //////Back substitution and iterative refiment + //////Back substitution and iterative refiment - //phase = 33; - //iparm[7] = 2; //max num of iterative refiment - //////set right hand b[i] - //PARDISO(pt, &maxfct, &mnum, &mtype,&phase, - // &n, spmat.get_values_array(), (int*) spmat.get_iofeir_array(), (int*) spmat.get_columns_array(), NULL, &nrhs, iparm, &msglvl, b, x, &error); - //if (error != 0) - //{ - // printf("\nERROR during solution: %d", error); - // exit(3); - //} - //printf("\nSolve completed..."); - //printf("\nThe solution of the system is:\n"); - ////output - //for (i=0;i