From 5ea84f25c94a63977fb3331d678a2a9543ce6656 Mon Sep 17 00:00:00 2001 From: akva Date: Thu, 8 Dec 2011 16:23:42 +0000 Subject: [PATCH] added: support for XML based input files git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1341 e10b68d5-8a6e-419e-a041-bce267b0401d --- src/LinAlg/LinSolParams.C | 49 +++++++++ src/LinAlg/LinSolParams.h | 6 ++ src/SIM/NonLinSIM.C | 61 +++++++++++ src/SIM/NonLinSIM.h | 5 + src/SIM/SIM1D.C | 158 ++++++++++++++++++++++++++++ src/SIM/SIM1D.h | 7 ++ src/SIM/SIM2D.C | 210 ++++++++++++++++++++++++++++++++++++++ src/SIM/SIM2D.h | 10 ++ src/SIM/SIM3D.C | 183 +++++++++++++++++++++++++++++++++ src/SIM/SIM3D.h | 12 +++ src/SIM/SIMbase.C | 204 ++++++++++++++++++++++++++++++++++++ src/SIM/SIMbase.h | 19 ++++ src/SIM/SIMinput.C | 66 +++++++++++- src/SIM/SIMinput.h | 16 +++ 14 files changed, 1005 insertions(+), 1 deletion(-) diff --git a/src/LinAlg/LinSolParams.C b/src/LinAlg/LinSolParams.C index 4220de9c..561a121f 100644 --- a/src/LinAlg/LinSolParams.C +++ b/src/LinAlg/LinSolParams.C @@ -18,6 +18,8 @@ #include "petscmg.h" #endif +#include "tinyxml.h" + void LinSolParams::setDefault () { @@ -147,6 +149,53 @@ bool LinSolParams::read (std::istream& is, int nparam) return true; } +static bool IsOK(const TiXmlElement* child, const char* value) +{ + return !strcasecmp(child->Value(),value) && + child->FirstChild() && child->FirstChild()->Value(); +} + +bool LinSolParams::read (const TiXmlElement* elem) +{ + const TiXmlElement *child = elem->FirstChildElement(); + while (child) { +#ifdef HAS_PETSC + if (IsOK(child,"type")) { + method = child->FirstChild()->Value(); + } else if (IsOK(child,"pc")) { + prec = child->FirstChild()->Value(); + } else if (IsOK(child,"hypretype")) { + hypretype = child->FirstChild()->Value(); + } else if (IsOK(child,"package")) { + package = child->FirstChild()->Value(); + } else if (IsOK(child,"levels")) { + levels = atoi(child->FirstChild()->Value()); + } else if (IsOK(child,"overlap")) { + overlap = atoi(child->FirstChild()->Value()); + } else if (IsOK(child,"atol")) { + atol = atof(child->FirstChild()->Value()); + } else if (IsOK(child,"rtol")) { + rtol = atof(child->FirstChild()->Value()); + } else if (IsOK(child,"dtol")) { + dtol = atof(child->FirstChild()->Value()); + } else if (IsOK(child,"maxits")) { + maxIts = atoi(child->FirstChild()->Value()); + } else if (IsOK(child,"nullspace")) { + if (!strcasecmp(child->FirstChild()->Value(),"constant")) + nullspc = CONSTANT; + else if (!strcasecmp(child->FirstChild()->Value(),"rigid_body")) + nullspc = RIGID_BODY; + } + else { + std::cerr <<" *** LinSolParams::read: Unknown keyword: " + << child->Value() << std::endl; + } +#endif + child = child->NextSiblingElement(); + } + return true; +} + bool LinSolParams::read (const char* filename) { diff --git a/src/LinAlg/LinSolParams.h b/src/LinAlg/LinSolParams.h index d704f7a6..fe091eac 100644 --- a/src/LinAlg/LinSolParams.h +++ b/src/LinAlg/LinSolParams.h @@ -23,6 +23,9 @@ #endif +class TiXmlElement; + + //! \brief Null-space for matrix operator enum NullSpace { NONE, CONSTANT, RIGID_BODY }; @@ -49,6 +52,9 @@ public: //! \brief Read linear solver parameters from stream. virtual bool read(std::istream& is, int nparams = 10); + //! \brief Read linear solver parameters from XML document. + virtual bool read(const TiXmlElement* elem); + //! \brief Read linear solver parameters from file. virtual bool read(const char* filename); diff --git a/src/SIM/NonLinSIM.C b/src/SIM/NonLinSIM.C index 527861fa..0b7b5c8d 100644 --- a/src/SIM/NonLinSIM.C +++ b/src/SIM/NonLinSIM.C @@ -17,6 +17,8 @@ #include "Utilities.h" #include +#include "tinyxml.h" + NonLinSIM::NonLinSIM (SIMbase* sim) : model(sim), nBlock(0) { @@ -111,6 +113,65 @@ void NonLinSIM::initSystem (SystemMatrix::Type mType, size_t nGauss) } +bool NonLinSIM::parse (const TiXmlElement* elem) +{ + if (strcasecmp(elem->Value(),"nonlinearsolver")) + return model->parse(elem); + + const TiXmlElement* child = elem->FirstChildElement(); + while (child) { + if (!strcasecmp(child->Value(),"timestepping")) { + const TiXmlElement* step = child->FirstChildElement("step"); + steps.clear(); + while (step) { + double dt; + std::pair,double> tstep; + double start, end; + if (step->Attribute("start")) + start = atof(step->Attribute("start")); + if (step->Attribute("end")) + end = atof(step->Attribute("end")); + if (steps.empty()) + startTime = start; + if (step->FirstChild() && step->FirstChild()->Value()) { + std::istringstream cline(child->FirstChild()->Value()); + cline >> dt; + if (dt > 1.0 && ceil(dt) == dt) { // number of steps specified + dt = (end-steps.empty()?start:steps.back().second)/dt; + tstep.first.push_back(dt); + } else while (!cline.fail() && !cline.bad()) { // steps specified + tstep.first.push_back(dt); + cline >> dt; + } + tstep.second = end; + steps.push_back(tstep); + } + } + stopTime = steps.back().second; + } else if (!strcasecmp(child->Value(),"maxits")) { + if (child->FirstChild() && child->FirstChild()->Value()) + maxit = atoi(child->FirstChild()->Value()); + } else if (!strcasecmp(child->Value(),"nupdate")) { + if (child->FirstChild() && child->FirstChild()->Value()) + nupdat = atoi(child->FirstChild()->Value()); + } else if (!strcasecmp(child->Value(),"rtol")) { + if (child->FirstChild() && child->FirstChild()->Value()) + convTol = atof(child->FirstChild()->Value()); + } else if (!strcasecmp(child->Value(),"dtol")) { + if (child->FirstChild() && child->FirstChild()->Value()) + divgLim = atof(child->FirstChild()->Value()); + } else if (!strcasecmp(child->Value(),"eta")) { + if (child->FirstChild() && child->FirstChild()->Value()) + eta = atof(child->FirstChild()->Value()); + } + + child = child->NextSiblingElement(); + } + + return true; +} + + void NonLinSIM::init (SolvePrm& param, const RealArray& initVal) { param.initTime(startTime,stopTime,steps); diff --git a/src/SIM/NonLinSIM.h b/src/SIM/NonLinSIM.h index bc05b5d4..6ecbb46a 100644 --- a/src/SIM/NonLinSIM.h +++ b/src/SIM/NonLinSIM.h @@ -171,6 +171,11 @@ public: //! \param is The file stream to read from virtual bool parse(char* keyWord, std::istream& is); + //! \brief Parses a data section from an input stream. + //! \param[in] keyWord Keyword of current data section to read + //! \param is The file stream to read from + virtual bool parse(const TiXmlElement* elem); + private: SIMbase* model; //!< The isogeometric FE model diff --git a/src/SIM/SIM1D.C b/src/SIM/SIM1D.C index eaec7c17..b83ca5fb 100644 --- a/src/SIM/SIM1D.C +++ b/src/SIM/SIM1D.C @@ -21,6 +21,8 @@ #include #include +#include "tinyxml.h" + SIM1D::SIM1D (unsigned char n_f) { @@ -218,6 +220,162 @@ bool SIM1D::parse (char* keyWord, std::istream& is) } +bool SIM1D::parseGeometryTag(const TiXmlElement* elem) +{ + // The remaining keywords are retained for backward compatibility with the + // prototype version. They enable direct specification of topology and + // properties as well as grid refinement without using the GPM module. + + if (!strcasecmp(elem->Value(),"refine")) { + ASMs1D* pch = 0; + bool uniform=true; + if (elem->Attribute("type") && !strcasecmp(elem->Attribute("type"),"nonuniform")) + uniform = false; + size_t lowpatch = 1, uppatch = 2; + + if (elem->Attribute("patch")) + lowpatch = uppatch = atoi(elem->Attribute("patch")); + if (elem->Attribute("lowerpatch")) { + lowpatch = atoi(elem->Attribute("lowerpatch")); + uppatch = myModel.size(); + } + if (elem->Attribute("upperpatch")) + uppatch = atoi(elem->Attribute("upperpatch")); + + if (lowpatch < 1 || uppatch > myModel.size()) + { + std::cerr <<" *** SIM1D::parse: Invalid patch index " + << "lower: " << lowpatch << " upper: " << uppatch << std::endl; + return false; + } + + if (uniform) + { + int addu=0; + if (elem->Attribute("u")) + addu = atoi(elem->Attribute("u")); + for (size_t j = lowpatch-1; j < uppatch; j++) { + if ((pch = dynamic_cast(myModel[j]))) { + std::cout <<"\tRefining P"<< j+1 + <<" "<< addu << std::endl; + pch->uniformRefine(addu); + } + } + } else { + char* cline2 = strdup(elem->FirstChildElement()->Value()); + char* cline = cline2; + strtok(cline," "); + RealArray xi; + while (cline) { + xi.push_back(atof(cline)); + cline = strtok(NULL," "); + } + free(cline2); + for (size_t j = lowpatch-1; j < uppatch; j++) { + if ((pch = dynamic_cast(myModel[j]))) + { + std::cout <<"\tRefining P"<< j+1; + for (size_t i = 0; i < xi.size(); i++) + std::cout <<" "<< xi[i]; + std::cout << std::endl; + pch->refine(xi); + } + } + } + } else if (!strcasecmp(elem->Value(),"raiseorder")) { + size_t lowpatch = 1, uppatch = 2; + if (elem->Attribute("patch")) + lowpatch = uppatch = atoi(elem->Attribute("patch")); + if (elem->Attribute("lowerpatch")) { + lowpatch = atoi(elem->Attribute("lowerpatch")); + uppatch = myModel.size(); + } + if (elem->Attribute("upperpatch")) + uppatch = atoi(elem->Attribute("upperpatch")); + + if (lowpatch < 1 || uppatch > myModel.size()) + { + std::cerr <<" *** SIM2D::parse: Invalid patch index " + << "lower: " << lowpatch << " upper: " << uppatch << std::endl; + return false; + } + int addu=0; + if (elem->Attribute("u")) + addu = atoi(elem->Attribute("u")); + ASMs1D* pch; + for (size_t j = lowpatch-1; j < uppatch; j++) { + if ((pch = dynamic_cast(myModel[j]))) { + std::cout <<"\tRaising order of P"<< j+1 <<" "<< addu << std::endl; + pch->raiseOrder(addu); + } + } + } + else if (!strcasecmp(elem->Value(),"topology")) { + if (createFEMmodel()) return false; + + const TiXmlElement* child = elem->FirstChildElement("connection"); + while (child) { + int master=0, slave=0, mVert=0, sVert=0; + if (child->Attribute("master")) + master = atoi(child->Attribute("master")); + if (child->Attribute("mvert")) + mVert = atoi(child->Attribute("mvert")); + if (child->Attribute("slave")) + slave = atoi(child->Attribute("slave")); + if (child->Attribute("svert")) + sVert = atoi(child->Attribute("svert")); + + if (master == slave || + master < 1 || master > (int)myModel.size() || + slave < 1 || slave > (int)myModel.size()) + { + std::cerr <<" *** SIM1D::parse: Invalid patch indices " + << master <<" "<< slave << std::endl; + return false; + } + std::cout <<"\tConnecting P"<< slave <<" V"<< sVert + <<" to P"<< master <<" E"<< mVert << std::endl; + ASMs1D* spch = static_cast(myModel[slave-1]); + ASMs1D* mpch = static_cast(myModel[master-1]); + if (!spch->connectPatch(sVert,*mpch,mVert)) + return false; + + child = child->NextSiblingElement(); + } + } else if (!strcasecmp(elem->Value(),"periodic")) { + if (createFEMmodel()) return false; + int patch=0; + if (elem->Attribute("patch")) + patch = atoi(elem->Attribute("patch")); + + if (patch < 1 || patch > (int)myModel.size()) { + std::cerr <<" *** SIM1D::parse: Invalid patch index " + << patch << std::endl; + return false; + } + std::cout <<"\tPeriodic P" << patch << std::endl; + static_cast(myModel[patch-1])->closeEnds(); + } + // TODO: constraints? fixpoints? + + return true; +} + + +bool SIM1D::parse(const TiXmlElement* elem) +{ + bool result=SIMbase::parse(elem); + + const TiXmlElement* child = elem->FirstChildElement(); + while (child) { + if (!strcasecmp(elem->Value(),"geometry")) + result &= parseGeometryTag(child); + child = child->NextSiblingElement(); + } + return result; +} + + /*! \brief Local-scope convenience function for error message generation. */ diff --git a/src/SIM/SIM1D.h b/src/SIM/SIM1D.h index b1109d44..e5ad2c71 100644 --- a/src/SIM/SIM1D.h +++ b/src/SIM/SIM1D.h @@ -48,6 +48,10 @@ protected: //! \param[in] pchInd 0-based index of the patch to read virtual bool readPatch(std::istream& isp, int pchInd); + //! \brief Parses a data section from an XML document. + //! \param[in] elem The XML element to parse + virtual bool parse(const TiXmlElement* elem); + //! \brief Preprocesses a user-defined Dirichlet boundary property. //! \param[in] patch 1-based index of the patch to receive the property //! \param[in] lndx Local index of the boundary item to receive the property @@ -59,6 +63,9 @@ protected: protected: unsigned char nf; //!< Number of scalar fields +private: + //! \brief Parse subtags of the block from the input file + bool parseGeometryTag(const TiXmlElement* elem); }; #endif diff --git a/src/SIM/SIM2D.C b/src/SIM/SIM2D.C index c49f0d99..bba896bd 100644 --- a/src/SIM/SIM2D.C +++ b/src/SIM/SIM2D.C @@ -19,6 +19,8 @@ #include #include +#include "tinyxml.h" + /*! A struct defining a patch interface for C1-continuous models. @@ -48,6 +50,179 @@ SIM2D::SIM2D (unsigned char n1, unsigned char n2) : isRefined(false) } +bool SIM2D::parseGeometryTag(const TiXmlElement* elem) +{ + // The remaining keywords are retained for backward compatibility with the + // prototype version. They enable direct specification of topology and + // properties as well as grid refinement without using the GPM module. + + if (!strcasecmp(elem->Value(),"refine") && !isRefined) { + ASM2D* pch = 0; + bool uniform=true; + if (elem->Attribute("type") && !strcasecmp(elem->Attribute("type"),"nonuniform")) + uniform = false; + size_t lowpatch = 1, uppatch = 2; + + if (elem->Attribute("patch")) + lowpatch = uppatch = atoi(elem->Attribute("patch")); + if (elem->Attribute("lowerpatch")) { + lowpatch = atoi(elem->Attribute("lowerpatch")); + uppatch = myModel.size(); + } + if (elem->Attribute("upperpatch")) + uppatch = atoi(elem->Attribute("upperpatch")); + + if (lowpatch < 1 || uppatch > myModel.size()) + { + std::cerr <<" *** SIM2D::parse: Invalid patch index " + << "lower: " << lowpatch << " upper: " << uppatch << std::endl; + return false; + } + + if (uniform) + { + int addu=0, addv = 0; + if (elem->Attribute("u")) + addu = atoi(elem->Attribute("u")); + if (elem->Attribute("v")) + addv = atoi(elem->Attribute("v")); + for (size_t j = lowpatch-1; j < uppatch; j++) { + if ((pch = dynamic_cast(myModel[j]))) { + std::cout <<"\tRefining P"<< j+1 + <<" "<< addu <<" "<< addv << std::endl; + pch->uniformRefine(0,addu); + pch->uniformRefine(1,addv); + } + } + } else { + int dir=1; + if (elem->Attribute("dir")) + dir = atoi(elem->Attribute("dir")); + char* cline2 = strdup(elem->FirstChildElement()->Value()); + char* cline = cline2; + strtok(cline," "); + RealArray xi; + while (cline) { + xi.push_back(atof(cline)); + cline = strtok(NULL," "); + } + free(cline2); + for (size_t j = lowpatch-1; j < uppatch; j++) { + if ((pch = dynamic_cast(myModel[j]))) + { + std::cout <<"\tRefining P"<< j+1 <<" dir="<< dir; + for (size_t i = 0; i < xi.size(); i++) + std::cout <<" "<< xi[i]; + std::cout << std::endl; + pch->refine(dir-1,xi); + } + } + } + } else if (!strcasecmp(elem->Value(),"raiseorder") && !isRefined) { + size_t lowpatch = 1, uppatch = 2; + if (elem->Attribute("patch")) + lowpatch = uppatch = atoi(elem->Attribute("patch")); + if (elem->Attribute("lowerpatch")) { + lowpatch = atoi(elem->Attribute("lowerpatch")); + uppatch = myModel.size(); + } + if (elem->Attribute("upperpatch")) + uppatch = atoi(elem->Attribute("upperpatch")); + + if (lowpatch < 1 || uppatch > myModel.size()) + { + std::cerr <<" *** SIM2D::parse: Invalid patch index " + << "lower: " << lowpatch << " upper: " << uppatch << std::endl; + return false; + } + int addu=0, addv = 0; + if (elem->Attribute("u")) + addu = atoi(elem->Attribute("u")); + if (elem->Attribute("v")) + addv = atoi(elem->Attribute("v")); + ASM2D* pch; + for (size_t j = lowpatch-1; j < uppatch; j++) { + if ((pch = dynamic_cast(myModel[j]))) { + std::cout <<"\tRaising order of P"<< j+1 + <<" "<< addu <<" "<< addv << std::endl; + pch->raiseOrder(addu,addv); + } + } + } + else if (!strcasecmp(elem->Value(),"topology")) { + if (createFEMmodel()) return false; + + const TiXmlElement* child = elem->FirstChildElement("connection"); + while (child) { + int master=0, slave=0, mEdge=0, sEdge=0; + if (child->Attribute("master")) + master = atoi(child->Attribute("master")); + if (child->Attribute("medge")) + mEdge = atoi(child->Attribute("medge")); + if (child->Attribute("slave")) + slave = atoi(child->Attribute("slave")); + if (child->Attribute("sedge")) + sEdge = atoi(child->Attribute("sedge")); + + bool rever=false; + if (child->Attribute("reverse") && !strcasecmp(child->Attribute("reverse"),"true")) + rever = true; + + if (master == slave || + master < 1 || master > (int)myModel.size() || + slave < 1 || slave > (int)myModel.size()) + { + std::cerr <<" *** SIM2D::parse: Invalid patch indices " + << master <<" "<< slave << std::endl; + return false; + } + std::cout <<"\tConnecting P"<< slave <<" E"<< sEdge + <<" to P"<< master <<" E"<< mEdge + <<" reversed? "<< rever << std::endl; + ASMs2D* spch = static_cast(myModel[slave-1]); + ASMs2D* mpch = static_cast(myModel[master-1]); + if (!spch->connectPatch(sEdge,*mpch,mEdge,rever)) + return false; + + child = child->NextSiblingElement(); + } + } else if (!strcasecmp(elem->Value(),"periodic")) { + if (createFEMmodel()) return false; + int patch=0, pedir=1; + if (elem->Attribute("patch")) + patch = atoi(elem->Attribute("patch")); + if (elem->Attribute("dir")) + pedir = atoi(elem->Attribute("dir")); + + if (patch < 1 || patch > (int)myModel.size()) { + std::cerr <<" *** SIM2D::parse: Invalid patch index " + << patch << std::endl; + return false; + } + std::cout <<"\tPeriodic "<< char('H'+pedir) <<"-direction P"<< patch + << std::endl; + static_cast(myModel[patch-1])->closeEdges(pedir); + } + // TODO: constraints? fixpoints? + + return true; +} + + +bool SIM2D::parse(const TiXmlElement* elem) +{ + bool result=SIMbase::parse(elem); + + const TiXmlElement* child = elem->FirstChildElement(); + while (child) { + if (!strcasecmp(elem->Value(),"geometry")) + result &= parseGeometryTag(child); + child = child->NextSiblingElement(); + } + return result; +} + + bool SIM2D::parse (char* keyWord, std::istream& is) { char* cline = 0; @@ -431,6 +606,41 @@ void SIM2D::clonePatches (const FEModelVec& patches, } + +void SIM2D::readNodes (std::istream& isn) +{ + while (isn.good()) { + int patch = 0; + isn >> patch; + ++patch; + int pid = getLocalPatchIndex(patch); + if (pid < 0) + return; + + ASMs2D::BlockNodes n; + for (size_t i = 0; i < 4 && isn.good(); i++) + isn >> n.ibnod[i]; + for (size_t i = 0; i < 4 && isn.good(); i++) + isn >> n.edges[i].icnod >> n.edges[i].incr; + isn >> n.iinod; + + // We always require the node numbers to be 1-based + for (size_t i = 0; i < 4; i++) + ++n.ibnod[i]; + for (size_t i = 0; i < 4; i++) + ++n.edges[i].icnod; + ++n.iinod; + + if (isn.good() && pid > 0) { + if (!static_cast(myModel[pid-1])->assignNodeNumbers(n)) { + std::cerr <<" *** SIM2D::readNodes: Failed to assign node numbers" + <<" for patch "<< patch << std::endl; + return; + } + } + } +} + bool SIM2D::refine (const std::vector& elements, const std::vector& options, const char* fName) { diff --git a/src/SIM/SIM2D.h b/src/SIM/SIM2D.h index 850462a6..a0442624 100644 --- a/src/SIM/SIM2D.h +++ b/src/SIM/SIM2D.h @@ -48,6 +48,10 @@ protected: //! \param is The file stream to read from virtual bool parse(char* keyWord, std::istream& is); + //! \brief Parses a data section from an XML document. + //! \param[in] elem The XML element to parse + virtual bool parse(const TiXmlElement* elem); + //! \brief Reads patches from given input stream. //! \param[in] isp The input stream to read from virtual bool readPatches(std::istream& isp); @@ -64,6 +68,10 @@ protected: virtual bool readNodes(std::istream& isn, int pchInd, int basis = 0, bool oneBased = false); + //! \brief Reads node numbers from given input stream. + //! \param[in] isn The file stream to read from + void readNodes(std::istream& isn); + //! \brief Refines a list of elements. //! \param[in] elements 1-based indices of the elements to refine //! \param[in] options Input options to refinement algorithm @@ -83,6 +91,8 @@ protected: private: bool isRefined; //!< If \e true, the model has been adaptively refined + //! \brief Parse subtags of the block from the input file + bool parseGeometryTag(const TiXmlElement* elem); protected: unsigned char nf[2]; //!< Number of scalar fields diff --git a/src/SIM/SIM3D.C b/src/SIM/SIM3D.C index 4be06351..11ac17ec 100644 --- a/src/SIM/SIM3D.C +++ b/src/SIM/SIM3D.C @@ -22,6 +22,8 @@ #include #include +#include "tinyxml.h" + SIM3D::SIM3D (bool checkRHS, unsigned char n1, unsigned char n2) { @@ -292,6 +294,148 @@ bool SIM3D::parse (char* keyWord, std::istream& is) } +bool SIM3D::parseGeometryTag(const TiXmlElement* elem) +{ + // The remaining keywords are retained for backward compatibility with the + // prototype version. They enable direct specification of topology and + // properties as well as grid refinement without using the GPM module. + + if (!strcasecmp(elem->Value(),"refine")) { + bool uniform=true; + if (elem->Attribute("type") && !strcasecmp(elem->Attribute("type"),"nonuniform")) + uniform = false; + size_t lowpatch = 1, uppatch = 2; + + if (elem->Attribute("patch")) + lowpatch = uppatch = atoi(elem->Attribute("patch")); + if (elem->Attribute("lowerpatch")) { + lowpatch = atoi(elem->Attribute("lowerpatch")); + uppatch = myModel.size(); + } + if (elem->Attribute("upperpatch")) + uppatch = atoi(elem->Attribute("upperpatch")); + + if (lowpatch < 1 || uppatch > myModel.size()) + { + std::cerr <<" *** SIM3D::parse: Invalid patch index " + << "lower: " << lowpatch << " upper: " << uppatch << std::endl; + return false; + } + + if (uniform) + { + int addu=0, addv = 0, addw = 0; + if (elem->Attribute("u")) + addu = atoi(elem->Attribute("u")); + if (elem->Attribute("v")) + addv = atoi(elem->Attribute("v")); + if (elem->Attribute("w")) + addw = atoi(elem->Attribute("w")); + for (size_t j = lowpatch-1; j < uppatch; j++) { + std::cout <<"\tRefining P"<< j+1 + <<" "<< addu <<" "<< addv <<" "<< addw << std::endl; + ASMs3D* pch = static_cast(myModel[j]); + pch->uniformRefine(0,addu); + pch->uniformRefine(1,addv); + pch->uniformRefine(2,addw); + } + } else { + int dir=1; + if (elem->Attribute("dir")) + dir = atoi(elem->Attribute("dir")); + char* cline2 = strdup(elem->FirstChildElement()->Value()); + char* cline = cline2; + strtok(cline," "); + RealArray xi; + while (cline) { + xi.push_back(atof(cline)); + cline = strtok(NULL," "); + } + free(cline2); + ASMs3D* pch; + for (size_t j = lowpatch-1; j < uppatch; j++) { + if ((pch = dynamic_cast(myModel[j]))) + { + std::cout <<"\tRefining P"<< j+1 <<" dir="<< dir; + for (size_t i = 0; i < xi.size(); i++) + std::cout <<" "<< xi[i]; + std::cout << std::endl; + pch->refine(dir-1,xi); + } + } + } + } else if (!strcasecmp(elem->Value(),"raiseorder")) { + size_t lowpatch = 1, uppatch = 2; + if (elem->Attribute("patch")) + lowpatch = uppatch = atoi(elem->Attribute("patch")); + if (elem->Attribute("lowerpatch")) { + lowpatch = atoi(elem->Attribute("lowerpatch")); + uppatch = myModel.size(); + } + if (elem->Attribute("upperpatch")) + uppatch = atoi(elem->Attribute("upperpatch")); + + if (lowpatch < 1 || uppatch > myModel.size()) + { + std::cerr <<" *** SIM3D::parse: Invalid patch index " + << "lower: " << lowpatch << " upper: " << uppatch << std::endl; + return false; + } + int addu=0, addv = 0, addw = 0; + if (elem->Attribute("u")) + addu = atoi(elem->Attribute("u")); + if (elem->Attribute("v")) + addv = atoi(elem->Attribute("v")); + if (elem->Attribute("w")) + addw = atoi(elem->Attribute("w")); + ASMs3D* pch; + for (size_t j = lowpatch-1; j < uppatch; j++) { + if ((pch = dynamic_cast(myModel[j]))) { + std::cout <<"\tRaising order of P"<< j+1 + <<" "<< addu <<" "<< addv <<" " << addw << std::endl; + pch->raiseOrder(addu,addv,addw); + } + } + } + else if (!strcasecmp(elem->Value(),"topology")) { // TODO? + assert(0); + } else if (!strcasecmp(elem->Value(),"periodic")) { + if (createFEMmodel()) return false; + int patch=0, pfdir=1; + if (elem->Attribute("patch")) + patch = atoi(elem->Attribute("patch")); + if (elem->Attribute("dir")) + pfdir = atoi(elem->Attribute("dir")); + + if (patch < 1 || patch > (int)myModel.size()) { + std::cerr <<" *** SIM3D::parse: Invalid patch index " + << patch << std::endl; + return false; + } + std::cout <<"\tPeriodic "<< char('H'+pfdir) <<"-direction P"<< patch + << std::endl; + static_cast(myModel[patch-1])->closeFaces(pfdir); + } + // TODO: topology? constraints? fixpoints? + + return true; +} + + +bool SIM3D::parse(const TiXmlElement* elem) +{ + bool result=SIMbase::parse(elem); + + const TiXmlElement* child = elem->FirstChildElement(); + while (child) { + if (!strcasecmp(elem->Value(),"geometry")) + result &= parseGeometryTag(child); + child = child->NextSiblingElement(); + } + return result; +} + + /*! \brief Local-scope convenience function for error message generation. */ @@ -528,3 +672,42 @@ bool SIM3D::readNodes (std::istream& isn, int pchInd, int basis, bool oneBased) return static_cast(myModel[pchInd])->assignNodeNumbers(n,basis); } + + +void SIM3D::readNodes (std::istream& isn) +{ + while (isn.good()) { + int patch = 0; + isn >> patch; + ++patch; + int pid = getLocalPatchIndex(patch); + if (pid < 0) + return; + + ASMs3D::BlockNodes n; + for (size_t i = 0; i < 8 && isn.good(); i++) + isn >> n.ibnod[i]; + for (size_t i = 0; i < 12 && isn.good(); i++) + isn >> n.edges[i].icnod >> n.edges[i].incr; + for (size_t i = 0; i < 6 && isn.good(); i++) + isn >> n.faces[i].isnod >> n.faces[i].incrI >> n.faces[i].incrJ; + isn >> n.iinod; + + // We always require the node numbers to be 1-based + for (size_t i = 0; i < 8; i++) + ++n.ibnod[i]; + for (size_t i = 0; i < 12; i++) + ++n.edges[i].icnod; + for (size_t i = 0; i < 6; i++) + ++n.faces[i].isnod; + ++n.iinod; + + if (isn.good() && pid > 0) { + if (!static_cast(myModel[pid-1])->assignNodeNumbers(n)) { + std::cerr <<" *** SIM3D::readNodes: Failed to assign node numbers" + <<" for patch "<< patch << std::endl; + return; + } + } + } +} diff --git a/src/SIM/SIM3D.h b/src/SIM/SIM3D.h index fd72178e..5296c22a 100644 --- a/src/SIM/SIM3D.h +++ b/src/SIM/SIM3D.h @@ -38,6 +38,10 @@ public: //! \param[in] ng Number of Gauss points in each parameter direction virtual void setQuadratureRule(size_t ng); + //! \brief Reads node numbers from given input stream. + //! \param[in] isn The file stream to read from + void readNodes(std::istream& isn); + protected: //! \brief Parses a data section from an input stream. //! \param[in] keyWord Keyword of current data section to read @@ -60,6 +64,11 @@ protected: virtual bool readNodes(std::istream& isn, int pchInd, int basis = 0, bool oneBased = false); + //! \brief Parses a data section from an XML element + //! \param[in] elem The XML element to parse + ////! \param is The file stream to read from + virtual bool parse(const TiXmlElement* elem); + //! \brief Preprocesses a user-defined Dirichlet boundary property. //! \param[in] patch 1-based index of the patch to receive the property //! \param[in] lndx Local index of the boundary item to receive the property @@ -81,6 +90,9 @@ protected: unsigned char nf[2]; //!< Number of scalar fields bool checkRHSys; //!< Check if all patches are in a right-hand system +private: + //! \brief Parse subtags of the block from the input file + bool parseGeometryTag(const TiXmlElement* elem); }; #endif diff --git a/src/SIM/SIMbase.C b/src/SIM/SIMbase.C index 9f3306a0..c7612c31 100644 --- a/src/SIM/SIMbase.C +++ b/src/SIM/SIMbase.C @@ -28,6 +28,7 @@ #include "Tensor.h" #include "Vec3Oper.h" #include "EigSolver.h" +#include "Functions.h" #include "Profiler.h" #include "ElementBlock.h" #include "VTF.h" @@ -38,6 +39,9 @@ #include #include #include +#include + +#include "tinyxml.h" ASM::Discretization SIMbase::discretization = ASM::Spline; @@ -114,6 +118,197 @@ void SIMbase::clearProperties () } +bool SIMbase::parseGeometryTag(const TiXmlElement* elem) +{ + if (!strcasecmp(elem->Value(),"patchfile")) { + if (myModel.empty()) { + std::string file = elem->FirstChild()->Value(); + std::cout <<"\nReading data file "<< file << std::endl; + std::ifstream isp(file.c_str()); + readPatches(isp); + + if (myModel.empty()) { + std::cerr <<" *** SIMbase::parse: No patches read"<< std::endl; + return false; + } + } + } else if (!strcasecmp(elem->Value(),"nodefile")) { + if (!createFEMmodel()) + return false; + std::string file = elem->FirstChild()->Value(); + std::ifstream isn(file.c_str()); + if (isn) { + std::cout <<"\nReading data file "<< file << std::endl; + } else { + std::cerr <<" *** SIMbase::read: Failure opening input file " + << file << std::endl; + return false; + } + readNodes(isn); + } else if (!strcasecmp(elem->Value(),"partitioning")) { + if (!elem->Attribute("procs") || atoi(elem->Attribute("procs")) != nProc) + return false; + if (myPid == 0) + std::cout <<"\nNumber of partitions: "<< nProc << std::endl; + + const TiXmlElement* part = elem->FirstChildElement("part"); + while (part) { + if (part->Attribute("proc") && atoi(part->Attribute("proc")) == myPid) { + int first=-2, last=-2; + if (part->Attribute("lower")) + first = atoi(part->Attribute("upper")); + if (part->Attribute("upper")) + last = atoi(part->Attribute("upper")); + + if (last > nGlPatches) nGlPatches = last; + + myPatches.reserve(last-first+1); + for (int j = first; j <= last && j > -1; j++) + myPatches.push_back(j); + } + part = part->NextSiblingElement("part"); + } + } + + return true; +} + + +bool SIMbase::parseOutputTag(const TiXmlElement* elem) +{ + if (!strcasecmp(elem->Value(),"resultpoints")) { + const TiXmlElement* point = elem->FirstChildElement("point"); + int i=0; + while (point) { + ResultPoint thePoint; + if (point->Attribute("patch")) + thePoint.patch = atoi(point->Attribute("patch")); + if (myPid == 0) + std::cout <<"\n\tPoint "<< i+1 <<": P"<< thePoint.patch <<" xi ="; + if (point->Attribute("u")) { + thePoint.par[0] = atof(point->Attribute("u")); + if (myPid == 0) + std::cout << ' ' << thePoint.par[0]; + } + if (point->Attribute("v")) { + thePoint.par[1] = atof(point->Attribute("v")); + if (myPid == 0) + std::cout << ' ' << thePoint.par[1]; + } + if (point->Attribute("w")) { + thePoint.par[2] = atoi(point->Attribute("w")); + if (myPid == 0) + std::cout << ' ' << thePoint.par[2]; + } + myPoints.push_back(thePoint); + i++; + point = point->NextSiblingElement(); + } + if (myPid == 0) + std::cout << std::endl; + } else if (!strcasecmp(elem->Value(),"vtfformat")) { + if (elem->FirstChild()->Value() && + !strcasecmp(elem->FirstChild()->Value(),"ascii")) + format = 0; + else if (elem->FirstChild()->Value() && + !strcasecmp(elem->FirstChild()->Value(),"binary")) + format = 1; + } else if (!strcasecmp(elem->Value(),"stride")) { + if (elem->FirstChild()->Value()) + vizIncr = atoi(elem->FirstChild()->Value()); + } + + return true; +} + + +bool SIMbase::parseBCTag(const TiXmlElement* elem) +{ + if (!strcasecmp(elem->Value(),"propertycodes")) { + const TiXmlElement* code = elem->FirstChildElement("code"); + while (code) { + int icode=0; + if (code->Attribute("value")) + icode = atoi(code->Attribute("value")); + const TiXmlElement* patch = code->FirstChildElement("patch"); + while (patch) { + Property p; + p.pindx = icode; + // TODO WTF is the only read number if lindex < 2 ? + if (patch->Attribute("face")) { + p.lindx = atoi(patch->Attribute("face")); + p.ldim = 2; + } + if (patch->Attribute("edge")) { + p.lindx = atoi(patch->Attribute("edge")); + p.ldim = 1; + } + if (patch->Attribute("vertex")) { + p.lindx = atoi(patch->Attribute("vertex")); + p.ldim = 0; + } + if (patch->Attribute("index")) + p.patch = getLocalPatchIndex(atoi(patch->Attribute("index"))); + if (p.patch > 0) + myProps.push_back(p); + patch = patch->NextSiblingElement("patch"); + } + code = code->NextSiblingElement(); + } + } else if (!strcasecmp(elem->Value(),"dirichlet")) { + if (ignoreDirichlet) // ignore all boundary conditions + return true; + int code=0; + if (elem->Attribute("code")) + code = atoi(elem->Attribute("code")); + double val=0.f; + if (elem->FirstChild() && elem->FirstChild()->Value()) + val = atof(elem->FirstChild()->Value()); + if (val == 0.f) { + setPropertyType(code,Property::DIRICHLET); + } else { + setPropertyType(code,Property::DIRICHLET_INHOM); + char* func=NULL; + if (elem->Attribute("function")) + func = strdup(elem->Attribute("function")); + if (func) { + char* func2 = func; + func = strtok(func," "); + myScalars[code] = const_cast(utl::parseRealFunc(const_cast(func),val)); + free(func2); + } else + myScalars[code] = new ConstFunc(val); + } + } + + return true; +} + + +bool SIMbase::parse (const TiXmlElement* elem) +{ + if (!strcasecmp(elem->Value(),"linearsolver")) + readLinSolParams(elem); + else { + const TiXmlElement* child = elem->FirstChildElement(); + bool result = true; + while (child) { + if (!strcasecmp(elem->Value(),"geometry")) + result &= parseGeometryTag(child); + else if (!strcasecmp(elem->Value(),"resultoutput")) + result &= parseOutputTag(child); + else if (!strcasecmp(elem->Value(),"boundaryconditions")) + result &= parseBCTag(child); + + child = child->NextSiblingElement(); + } + return result; + } + + return true; +} + + bool SIMbase::parse (char* keyWord, std::istream& is) { char* cline = 0; @@ -358,6 +553,15 @@ void SIMbase::readLinSolParams (std::istream& is, int npar) } +void SIMbase::readLinSolParams (const TiXmlElement* elem) +{ + if (!mySolParams) + mySolParams = new LinSolParams(); + + mySolParams->read(elem); +} + + bool SIMbase::createFEMmodel () { for (size_t i = 0; i < myModel.size(); i++) diff --git a/src/SIM/SIMbase.h b/src/SIM/SIMbase.h index ea391244..90ad6369 100644 --- a/src/SIM/SIMbase.h +++ b/src/SIM/SIMbase.h @@ -97,6 +97,10 @@ public: //! \param is The file stream to read from virtual bool parse(char* keyWord, std::istream& is); + //! \brief Parses a data section from an xml document. + //! \param[in] elem The XML element to parse + virtual bool parse(const TiXmlElement* elem); + //! \brief Refines a list of elements. virtual bool refine(const std::vector&, const std::vector&, const char* = 0) { return false; } @@ -484,6 +488,9 @@ protected: //! \brief Creates the computational FEM model from the spline patches. bool createFEMmodel(); + //! \brief Reads node numbers from given input stream. + //! \param[in] isp The file stream to read from + virtual void readNodes(std::istream& isp) {} public: //! \brief Returns the local patch index for the given global patch number. //! \details For serial applications this is an identity mapping only, whereas @@ -554,6 +561,10 @@ protected: //! \details This method helps with encapsulating PETSc in libIFEM. void readLinSolParams(std::istream& is, int npar); + //! \brief Reads a LinSolParams object from the given XML element + //! \details This method helps with encapsulating PETSc in libIFEM. + void readLinSolParams(const TiXmlElement* elem); + //! \brief Finalizes the global equation system assembly. virtual bool finalizeAssembly(bool newLHSmatrix = true); @@ -603,6 +614,14 @@ protected: AlgEqSystem* myEqSys; //!< The actual linear equation system SAMpatch* mySam; //!< Auxiliary data for FE assembly management LinSolParams* mySolParams; //!< Input parameters for PETSc +private: + //! \brief Parses a subelement of the tag from input file + //! \param[in] elem The XML element to parse + bool parseGeometryTag(const TiXmlElement* elem); + //! \brief Parses a subelement of the tag from input file + //! \param[in] elem The XML element to parse + bool parseOutputTag(const TiXmlElement* elem); + bool parseBCTag(const TiXmlElement* elem); }; #endif diff --git a/src/SIM/SIMinput.C b/src/SIM/SIMinput.C index 96f03c11..2109dc26 100644 --- a/src/SIM/SIMinput.C +++ b/src/SIM/SIMinput.C @@ -18,6 +18,8 @@ #endif #include +#include "tinyxml.h" + int SIMinput::msgLevel = 2; @@ -34,7 +36,16 @@ SIMinput::SIMinput () } -bool SIMinput::read (const char* fileName) +bool SIMinput::read(const char* fileName) +{ + if (strcasestr(fileName,".xinp")) + return readXML(fileName); + + return readFlat(fileName); +} + + +bool SIMinput::readFlat (const char* fileName) { std::ifstream is(fileName); if (is) @@ -58,3 +69,56 @@ bool SIMinput::read (const char* fileName) std::cout <<"\nReading input file succeeded."<< std::endl; return true; } + + +bool SIMinput::readXML (const char* fileName) +{ + TiXmlDocument doc; + if (!doc.LoadFile(fileName)) { + std::cerr << __PRETTY_FUNCTION__ << ": Failed to load " << fileName << std::endl; + std::cerr << "\t Error at line " << doc.ErrorRow() << ": " << doc.ErrorDesc() << std::endl; + return false; + } + + if (!doc.RootElement() || + strcmp(doc.RootElement()->Value(),"simulation")) { + std::cerr << __PRETTY_FUNCTION__ << ": Malformatted input file " << fileName << std::endl; + return false; + } + + std::vector parsed = handlePriorityTags(doc.RootElement()); + // now parse the rest + TiXmlElement* elem = doc.RootElement()->FirstChildElement(); + while (elem) { + if (find(parsed.begin(),parsed.end(),elem) == parsed.end()) { + if (!parse(elem)) { + std::cerr <<" *** SIMinput::read: Failure occured while parsing \"" + << elem->Value() <<"\""<< std::endl; + return false; + } + } + elem = elem->NextSiblingElement(); + } + return true; +} + +std::vector SIMinput::handlePriorityTags(const TiXmlElement* base) +{ + // add particular keywords to this list. these will be parsed first, + // in the order specified in the list + const char* special[] = {"geometry","boundaryconditions"}; + + std::vector parsed; + for (size_t i=0;iFirstChildElement(special[i]); + if (elem) { + if (!parse(elem)) { + std::cerr <<" *** SIMinput::read: Failure occured while parsing \"" + << elem->Value() <<"\""<< std::endl; + } + parsed.push_back(elem); + } + } + + return parsed; +} diff --git a/src/SIM/SIMinput.h b/src/SIM/SIMinput.h index c721c52d..aab4d243 100644 --- a/src/SIM/SIMinput.h +++ b/src/SIM/SIMinput.h @@ -15,12 +15,15 @@ #define _SIM_INPUT_H #include +#include /*! \brief Base class for NURBS-based FEM simulators with input file parsing. */ +class TiXmlElement; + class SIMinput { protected: @@ -33,14 +36,27 @@ public: //! \brief Reads model data from the specified input file \a *fileName. virtual bool read(const char* fileName); + //! \brief Reads a flat text input file + virtual bool readFlat(const char* fileName); + //! \brief Reads an XML input file + bool readXML(const char* fileName); //! \brief Parses a data section from an input stream. virtual bool parse(char* keyWord, std::istream& is) = 0; + //! \brief Parses a data section from an XML element. + virtual bool parse(const TiXmlElement* elem) { return false; } static int msgLevel; //!< Controls the amount of console output during solving protected: int myPid; //!< Processor ID in parallel simulations int nProc; //!< Number of processors in parallel simulations + + //! \brief Certain tags needs to be parsed before others. This function takes care of this + //! \param[in] base The base tag containing the elements to be prioritized + //! \returns A vector with elements that was handled + //! \details This function needs to be called in the app specific SIM classes + // prior parsing its app specific block. + std::vector handlePriorityTags(const TiXmlElement* base); }; #endif