From d71c9c1c50577ad48ecc5f0fcd1ea215674f71d3 Mon Sep 17 00:00:00 2001 From: kmo Date: Tue, 13 Mar 2012 14:41:55 +0000 Subject: [PATCH] Moved parsing of expression functions for analytic solutions to two new AnaSol constructors git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1520 e10b68d5-8a6e-419e-a041-bce267b0401d --- Apps/LinearElasticity/SIMLinEl2D.C | 61 ++------- Apps/LinearElasticity/SIMLinEl3D.C | 61 ++------- Apps/LinearElasticity/SIMLinElBeamC1.C | 137 ++++++++++--------- Apps/LinearElasticity/SIMLinElKL.C | 174 ++++++++++++------------- src/Integrands/AnaSol.C | 108 +++++++++++++++ src/Integrands/AnaSol.h | 13 ++ 6 files changed, 302 insertions(+), 252 deletions(-) create mode 100644 src/Integrands/AnaSol.C diff --git a/Apps/LinearElasticity/SIMLinEl2D.C b/Apps/LinearElasticity/SIMLinEl2D.C index ceb015b1..8a192d61 100644 --- a/Apps/LinearElasticity/SIMLinEl2D.C +++ b/Apps/LinearElasticity/SIMLinEl2D.C @@ -19,8 +19,6 @@ #include "Utilities.h" #include "Property.h" #include "AnaSol.h" -#include - #include "tinyxml.h" @@ -106,8 +104,8 @@ bool SIMLinEl2D::parse (char* keyWord, std::istream& is) else if (!strncasecmp(keyWord,"ANASOL",6)) { - cline = strtok(keyWord+6," "); int code = -1; + cline = strtok(keyWord+6," "); if (!strncasecmp(cline,"HOLE",4)) { double a = atof(strtok(NULL," ")); @@ -155,32 +153,10 @@ bool SIMLinEl2D::parse (char* keyWord, std::istream& is) } else if (!strncasecmp(cline,"EXPRESSION",10)) { - std::string variables, stress; - int lines = atoi(strtok(NULL, " ")); - char* c = strtok(NULL, " "); - if (c) - code = atoi(c); - else - code = 0; - for (int i = 0; i < lines; i++) { - std::string function = utl::readLine(is); - size_t pos = function.find("Variables="); - if (pos != std::string::npos) { - variables += function.substr(pos+10); - if (variables[variables.size()-1] != ';') - variables += ";"; - } - pos = function.find("Stress="); - if (pos != std::string::npos) { - stress = function.substr(pos+7); - mySol = new AnaSol(new STensorFuncExpr(stress,variables)); - std::cout <<"\nAnalytical solution:"; - if (!variables.empty()) - std::cout <<"\n\tVariables = "<< variables; - std::cout <<"\n\tStress = "<< stress << std::endl; - break; - } - } + std::cout <<"\nAnalytical solution: Expression"<< std::endl; + int lines = (cline = strtok(NULL," ")) ? atoi(cline) : 0; + code = (cline = strtok(NULL," ")) ? atoi(cline) : 0; + mySol = new AnaSol(is,lines,false); } else { @@ -339,9 +315,9 @@ bool SIMLinEl2D::parse (const TiXmlElement* elem) !strcasecmp(child->Value(),"linearpressure")) { if (child->FirstChild() && child->FirstChild()->Value()) { double p = atof(child->FirstChild()->Value()); - int code = 0, pdir = 1; - utl::getAttribute(child,"code",code); - utl::getAttribute(child,"dir",pdir); + int code = 0, pdir = 1; + utl::getAttribute(child,"code",code); + utl::getAttribute(child,"dir",pdir); setPropertyType(code,Property::NEUMANN); if (!strcasecmp(child->Value(),"linearpressure")) { RealFunc* pfl = new ConstTimeFunc(new LinearFunc(p)); @@ -349,7 +325,7 @@ bool SIMLinEl2D::parse (const TiXmlElement* elem) } else myTracs[code] = new PressureField(p,pdir); - if (myPid == 0) + if (myPid == 0) std::cout <<"\tPressure code "<< code <<" direction "<< pdir <<": "<< p << std::endl; } @@ -405,25 +381,12 @@ bool SIMLinEl2D::parse (const TiXmlElement* elem) <<" u0="<< u0 <<" E="<< E << std::endl; } else if (type == "expression") { - std::string variables, stress; - const TiXmlElement* var = child->FirstChildElement("variables"); - if (var && var->FirstChild() && var->FirstChild()->Value()) { - variables = var->FirstChild()->Value(); - if (variables[variables.size()-1] != ';') - variables += ";"; - } - const TiXmlElement* str = child->FirstChildElement("stress"); - if (str && str->FirstChild() && str->FirstChild()->Value()) - stress = str->FirstChild()->Value(); - mySol = new AnaSol(new STensorFuncExpr(stress,variables)); - std::cout <<"\nAnalytical solution:"; - if (!variables.empty()) - std::cout <<"\n\tVariables = "<< variables; - std::cout <<"\n\tStress = "<< stress << std::endl; + std::cout <<"\nAnalytical solution: Expression"<< std::endl; + mySol = new AnaSol(child,false); } else std::cerr <<" ** SIMLinEl2D::parse: Unknown analytical solution " - << type << std::endl; + << type <<" (ignored)"<< std::endl; // Define the analytical boundary traction field int code = 0; diff --git a/Apps/LinearElasticity/SIMLinEl3D.C b/Apps/LinearElasticity/SIMLinEl3D.C index 86a25ea7..c71ace77 100644 --- a/Apps/LinearElasticity/SIMLinEl3D.C +++ b/Apps/LinearElasticity/SIMLinEl3D.C @@ -219,8 +219,8 @@ bool SIMLinEl3D::parse (char* keyWord, std::istream& is) else if (!strncasecmp(keyWord,"ANASOL",6)) { - cline = strtok(keyWord+6," "); int code = -1; + cline = strtok(keyWord+6," "); if (!strncasecmp(cline,"HOLE",4)) { double a = atof(strtok(NULL," ")); @@ -250,32 +250,10 @@ bool SIMLinEl3D::parse (char* keyWord, std::istream& is) } else if (!strncasecmp(cline,"EXPRESSION",10)) { - std::string variables, stress; - int lines = atoi(strtok(NULL, " ")); - char* c = strtok(NULL, " "); - if (c) - code = atoi(c); - else - code = 0; - for (int i = 0; i < lines; i++) { - std::string function = utl::readLine(is); - size_t pos = function.find("Variables="); - if (pos != std::string::npos) { - variables += function.substr(pos+10); - if (variables[variables.size()-1] != ';') - variables += ";"; - } - pos = function.find("Stress="); - if (pos != std::string::npos) { - stress = function.substr(pos+7); - mySol = new AnaSol(new STensorFuncExpr(stress,variables)); - std::cout <<"\nAnalytical solution:"; - if (!variables.empty()) - std::cout <<"\n\tVariables = "<< variables; - std::cout <<"\n\tStress = "<< stress << std::endl; - break; - } - } + std::cout <<"\nAnalytical solution: Expression"<< std::endl; + int lines = (cline = strtok(NULL," ")) ? atoi(cline) : 0; + code = (cline = strtok(NULL," ")) ? atoi(cline) : 0; + mySol = new AnaSol(is,lines,false); } else { @@ -475,9 +453,9 @@ bool SIMLinEl3D::parse (const TiXmlElement* elem) !strcasecmp(child->Value(),"linearpressure")) { if (child->FirstChild() && child->FirstChild()->Value()) { double p = atof(child->FirstChild()->Value()); - int code = 0, pdir = 1; - utl::getAttribute(child,"code",code); - utl::getAttribute(child,"dir",pdir); + int code = 0, pdir = 1; + utl::getAttribute(child,"code",code); + utl::getAttribute(child,"dir",pdir); setPropertyType(code,Property::NEUMANN); if (!strcasecmp(child->Value(),"linearpressure")) { RealFunc* pfl = new ConstTimeFunc(new LinearFunc(p)); @@ -485,7 +463,7 @@ bool SIMLinEl3D::parse (const TiXmlElement* elem) } else myTracs[code] = new PressureField(p,pdir); - if (myPid == 0) + if (myPid == 0) std::cout <<"\tPressure code "<< code <<" direction "<< pdir <<": "<< p << std::endl; } @@ -522,30 +500,17 @@ bool SIMLinEl3D::parse (const TiXmlElement* elem) <<" F0="<< F0 << std::endl; } else if (type == "expression") { - std::string variables, stress; - const TiXmlElement* var = child->FirstChildElement("variables"); - if (var && var->FirstChild() && var->FirstChild()->Value()) { - variables = var->FirstChild()->Value(); - if (variables[variables.size()-1] != ';') - variables += ";"; - } - const TiXmlElement* str = child->FirstChildElement("stress"); - if (str && str->FirstChild() && str->FirstChild()->Value()) - stress = str->FirstChild()->Value(); - mySol = new AnaSol(new STensorFuncExpr(stress,variables)); - std::cout <<"\nAnalytical solution:"; - if (!variables.empty()) - std::cout <<"\n\tVariables = "<< variables; - std::cout <<"\n\tStress = "<< stress << std::endl; + std::cout <<"\nAnalytical solution: Expression"<< std::endl; + mySol = new AnaSol(child,false); } else std::cerr <<" ** SIMLinEl3D::parse: Unknown analytical solution " - << type << std::endl; + << type <<" (ignored)"<< std::endl; // Define the analytical boundary traction field int code = 0; utl::getAttribute(child,"code",code); - if (code > 0 && mySol &&mySol->getStressSol()) { + if (code > 0 && mySol && mySol->getStressSol()) { std::cout <<"Pressure code "<< code <<": Analytical traction"<< std::endl; setPropertyType(code,Property::NEUMANN); diff --git a/Apps/LinearElasticity/SIMLinElBeamC1.C b/Apps/LinearElasticity/SIMLinElBeamC1.C index e6cf87fa..3e41c11f 100644 --- a/Apps/LinearElasticity/SIMLinElBeamC1.C +++ b/Apps/LinearElasticity/SIMLinElBeamC1.C @@ -18,13 +18,11 @@ #include "AlgEqSystem.h" #include "ASMbase.h" #include "SAMpatch.h" -#include "Utilities.h" #include "Functions.h" +#include "Utilities.h" #include "Property.h" #include "AnaSol.h" #include "Vec3Oper.h" -#include - #include "tinyxml.h" @@ -112,12 +110,13 @@ bool SIMLinElBeamC1::parse (char* keyWord, std::istream& is) this->setPropertyType(code,Property::BODYLOAD); } } - /* + else if (!strncasecmp(keyWord,"ANASOL",6)) { if (mySol) return true; cline = strtok(keyWord+6," "); + /* if (!strncasecmp(cline,"NAVIERBEAM",7)) { double a = atof(strtok(NULL," ")); @@ -142,14 +141,21 @@ bool SIMLinElBeamC1::parse (char* keyWord, std::istream& is) else mySol = new AnaSol(new NavierBeam(a,t,E,pz)); } + */ + if (!strncasecmp(cline,"EXPRESSION",10)) + { + std::cout <<"\nAnalytical solution: Expression"<< std::endl; + int lines = (cline = strtok(NULL," ")) ? atoi(cline) : 0; + mySol = new AnaSol(is,lines,false); + } else { std::cerr <<" ** SIMLinElBeamC1::parse: Unknown analytical solution " - << cline <<" (ignored)"<< std::endl; + << cline <<" (ignored)"<< std::endl; return true; } } - */ + else return this->SIM1D::parse(keyWord,is); @@ -160,85 +166,88 @@ bool SIMLinElBeamC1::parse (char* keyWord, std::istream& is) bool SIMLinElBeamC1::parse(const TiXmlElement* elem) { if (strcasecmp(elem->Value(),"eulerbernoulli")) - return SIM1D::parse(elem); + return this->SIM1D::parse(elem); KirchhoffLovePlate* klp = dynamic_cast(myProblem); - if (!klp) - return false; + if (!klp) return false; std::vector parsed = handlePriorityTags(elem); + const TiXmlElement* child = elem->FirstChildElement(); while (child) { - if (find(parsed.begin(),parsed.end(),child) != parsed.end()) { - child = child->NextSiblingElement(); - continue; - } - if (!strcasecmp(child->Value(),"gravity")) { - double g=0; - if (child->Attribute("g")) - g = atof(child->Attribute("g")); - klp->setGravity(g); + if (find(parsed.begin(),parsed.end(),child) != parsed.end()) + ; + + else if (!strcasecmp(child->Value(),"gravity")) { + double g = 0.0; + utl::getAttribute(child,"g",g); if (myPid == 0) std::cout <<"\nGravitation constant: "<< g << std::endl; - } else if (!strcasecmp(child->Value(),"isotropic")) { - int code=0; - if (child->Attribute("code")) - code = atoi(child->Attribute("code")); - if (code > 0) + klp->setGravity(g); + } + + else if (!strcasecmp(child->Value(),"isotropic")) { + int code = 0; + if (utl::getAttribute(child,"code",code) && code > 0) setPropertyType(code,Property::MATERIAL,mVec.size()); - double E=1000.f, rho = 1.f, thk=0; - if (child->Attribute("E")) - E = atof(child->Attribute("E")); - if (child->Attribute("rho")) - rho = atof(child->Attribute("rho")); - if (child->Attribute("thickness")) - thk = atof(child->Attribute("thickness")); + double E = 1000.0, rho = 1.0, thk = 0.1; + utl::getAttribute(child,"E",E); + utl::getAttribute(child,"rho",rho); + utl::getAttribute(child,"thickness",thk); mVec.push_back(new LinIsotropic(E,0.0,rho,true)); tVec.push_back(thk); if (myPid == 0) std::cout <<"\tMaterial code "<< code <<": " - << E <<" " << rho <<" " << thk << std::endl; - if (!mVec.empty()) - klp->setMaterial(mVec.front()); - if (!tVec.empty() && tVec.front() != 0.0) + << E <<" "<< rho <<" " << thk << std::endl; + klp->setMaterial(mVec.front()); + if (tVec.front() != 0.0) klp->setThickness(tVec.front()); - } else if (!strcasecmp(child->Value(),"pointload")) { - PointLoad load; - if (child->Attribute("patch")) - load.patch = atoi(child->Attribute("patch")); - if (child->Attribute("xi")) - load.xi = atof(child->Attribute("xi")); - if (child->FirstChild() && child->FirstChild()->Value()) + } + + else if (!strcasecmp(child->Value(),"pointload")) { + PointLoad load; int patch; + utl::getAttribute(child,"patch",patch); + if (child->FirstChild()) { + load.patch = patch; load.pload = atof(child->FirstChild()->Value()); - if (myPid == 0) - std::cout <<"\n\tPoint: P" << load.patch - <<" xi = " << load.xi - <<" load = " << load.pload; - myLoads.push_back(load); - } else if (!strcasecmp(child->Value(),"pressure")) { - int code=0; - double p=0; - if (child->Attribute("code")) - code = atoi(child->Attribute("code")); - if (child->FirstChild() && child->FirstChild()->Value()) + utl::getAttribute(child,"xi",load.xi); + if (myPid == 0) + std::cout <<"\n\tPoint: P"<< load.patch + <<" xi = "<< load.xi + <<" load = "<< load.pload; + myLoads.push_back(load); + } + } + + else if (!strcasecmp(child->Value(),"pressure")) { + int code = 0; + double p = 0.0; + char* fn = NULL; + utl::getAttribute(child,"code",code); + if (child->FirstChild()) p = atof(child->FirstChild()->Value()); if (myPid == 0) std::cout <<"\tPressure code "<< code <<": "; - char* function; - if (child->Attribute("function")) - function = strdup(child->Attribute("function")); - else - function = strdup(""); - char* function2 = function; - function = strtok(function," "); - myScalars[code] = const_cast(utl::parseRealFunc(const_cast(function),p)); - free(function2); + std::string function; + if (utl::getAttribute(child,"function",function)) + fn = strtok(const_cast(function.c_str())," "); + myScalars[code] = const_cast(utl::parseRealFunc(fn,p)); std::cout << std::endl; if (code > 0) - setPropertyType(code,Property::BODYLOAD); + setPropertyType(code,Property::BODYLOAD); + } + + else if (!strcasecmp(child->Value(),"anasol")) { + std::string type; + utl::getAttribute(child,"type",type,true); + if (type == "expression") { + std::cout <<"\nAnalytical solution: Expression"<< std::endl; + mySol = new AnaSol(child); + } + else + std::cerr <<" ** SIMLinElKL::parse: Unknown analytical solution " + << type <<" (ignored)"<< std::endl; } - else - SIM1D::parse(child); child = child->NextSiblingElement(); } diff --git a/Apps/LinearElasticity/SIMLinElKL.C b/Apps/LinearElasticity/SIMLinElKL.C index 6f6aae60..914326d3 100644 --- a/Apps/LinearElasticity/SIMLinElKL.C +++ b/Apps/LinearElasticity/SIMLinElKL.C @@ -18,13 +18,11 @@ #include "AlgEqSystem.h" #include "ASMbase.h" #include "SAMpatch.h" -#include "Utilities.h" #include "Functions.h" +#include "Utilities.h" #include "Property.h" #include "AnaSol.h" #include "Vec3Oper.h" -#include - #include "tinyxml.h" @@ -149,6 +147,12 @@ bool SIMLinElKL::parse (char* keyWord, std::istream& is) else mySol = new AnaSol(new NavierPlate(a,b,t,E,nu,pz)); } + else if (!strncasecmp(cline,"EXPRESSION",10)) + { + std::cout <<"\nAnalytical solution: Expression"<< std::endl; + int lines = (cline = strtok(NULL," ")) ? atoi(cline) : 0; + mySol = new AnaSol(is,lines,false); + } else { std::cerr <<" ** SIMLinElKL::parse: Unknown analytical solution " @@ -164,132 +168,120 @@ bool SIMLinElKL::parse (char* keyWord, std::istream& is) } -bool SIMLinElKL::parse(const TiXmlElement* elem) +bool SIMLinElKL::parse (const TiXmlElement* elem) { if (strcasecmp(elem->Value(),"kirchhofflove")) - return SIMLinEl2D::parse(elem); + return this->SIMLinEl2D::parse(elem); KirchhoffLovePlate* klp = dynamic_cast(myProblem); - if (!klp) - return false; + if (!klp) return false; std::vector parsed = handlePriorityTags(elem); + const TiXmlElement* child = elem->FirstChildElement(); while (child) { - if (find(parsed.begin(),parsed.end(),child) != parsed.end()) { - child = child->NextSiblingElement(); - continue; - } - if (!strcasecmp(child->Value(),"gravity")) { - double g=0; - if (child->Attribute("g")) - g = atof(child->Attribute("g")); + if (find(parsed.begin(),parsed.end(),child) != parsed.end()) + ; + + else if (!strcasecmp(child->Value(),"gravity")) { + double g = 0.0; + utl::getAttribute(child,"g",g); klp->setGravity(g); if (myPid == 0) std::cout <<"\nGravitation constant: "<< g << std::endl; - } else if (!strcasecmp(child->Value(),"isotropic")) { - int code=0; - if (child->Attribute("code")) - code = atoi(child->Attribute("code")); - if (code > 0) + } + + else if (!strcasecmp(child->Value(),"isotropic")) { + int code = 0; + if (utl::getAttribute(child,"code",code) && code > 0) setPropertyType(code,Property::MATERIAL,mVec.size()); - double E=1000.f, nu=0.3f, rho = 1.f, thk=0; - if (child->Attribute("E")) - E = atof(child->Attribute("E")); - if (child->Attribute("rho")) - rho = atof(child->Attribute("rho")); - if (child->Attribute("nu")) - nu = atof(child->Attribute("nu")); - if (child->Attribute("thickness")) - thk = atof(child->Attribute("thickness")); + double E = 1000.0, nu = 0.3, rho = 1.0, thk = 0.1; + utl::getAttribute(child,"E",E); + utl::getAttribute(child,"nu",nu); + utl::getAttribute(child,"rho",rho); + utl::getAttribute(child,"thickness",thk); mVec.push_back(new LinIsotropic(E,nu,rho,true)); tVec.push_back(thk); if (myPid == 0) std::cout <<"\tMaterial code "<< code <<": " << E <<" "<< nu <<" "<< rho <<" " << thk << std::endl; - if (!mVec.empty()) - klp->setMaterial(mVec.front()); - if (!tVec.empty() && tVec.front() != 0.0) + klp->setMaterial(mVec.front()); + if (tVec.front() != 0.0) klp->setThickness(tVec.front()); - } else if (!strcasecmp(child->Value(),"pointload")) { - PointLoad load; - if (child->Attribute("patch")) - load.patch = atoi(child->Attribute("patch")); - if (child->Attribute("xi")) - load.xi[0] = atof(child->Attribute("xi")); - if (child->Attribute("eta")) - load.xi[1] = atof(child->Attribute("eta")); - if (child->FirstChild() && child->FirstChild()->Value()) + } + + else if (!strcasecmp(child->Value(),"pointload")) { + PointLoad load; int patch; + utl::getAttribute(child,"patch",patch); + if (child->FirstChild()) { + load.patch = patch; load.pload = atof(child->FirstChild()->Value()); + utl::getAttribute(child,"xi",load.xi[0]); + utl::getAttribute(child,"eta",load.xi[1]); if (myPid == 0) std::cout <<"\n\tPoint: P"<< load.patch <<" xi = "<< load.xi[0] <<" "<< load.xi[1] <<" load = "<< load.pload; myLoads.push_back(load); - } else if (!strcasecmp(child->Value(),"pressure")) { - int code=0; - double p=0; - if (child->Attribute("code")) - code = atoi(child->Attribute("code")); - if (child->FirstChild() && child->FirstChild()->Value()) + } + } + + else if (!strcasecmp(child->Value(),"pressure")) { + int code = 0; + double p = 0.0; + char* fn = NULL; + utl::getAttribute(child,"code",code); + if (child->FirstChild()) p = atof(child->FirstChild()->Value()); if (myPid == 0) std::cout <<"\tPressure code "<< code <<": "; - char* function; - if (child->Attribute("function")) - function = strdup(child->Attribute("function")); - else - function = strdup(""); - char* function2 = function; - function = strtok(function," "); - myScalars[code] = const_cast(utl::parseRealFunc(const_cast(function),p)); - free(function2); + std::string function; + if (utl::getAttribute(child,"function",function)) + fn = strtok(const_cast(function.c_str())," "); + myScalars[code] = const_cast(utl::parseRealFunc(fn,p)); std::cout << std::endl; if (code > 0) - setPropertyType(code,Property::BODYLOAD); - } else if (!strcasecmp(child->Value(),"anasol")) { - if (child->Attribute("type") && - !strcasecmp(child->Attribute("type"),"navierplate")) { - double a=0, b=0, c=0, d=0, t=0, E=0, nu=0, pz=0, xi=0, eta=0; - if (child->Attribute("a")) - a = atof(child->Attribute("a")); - if (child->Attribute("b")) - b = atof(child->Attribute("b")); - if (child->Attribute("c")) - c = atof(child->Attribute("c")); - if (child->Attribute("d")) - d = atof(child->Attribute("d")); - if (child->Attribute("t")) - t = atof(child->Attribute("t")); - if (child->Attribute("E")) - E = atof(child->Attribute("E")); - if (child->Attribute("nu")) - nu = atof(child->Attribute("nu")); - if (child->Attribute("pz")) - pz = atof(child->Attribute("pz")); - if (child->Attribute("xi")) - xi = atof(child->Attribute("xi")); - if (child->Attribute("eta")) - eta = atof(child->Attribute("eta")); + setPropertyType(code,Property::BODYLOAD); + } + + else if (!strcasecmp(child->Value(),"anasol")) { + std::string type; + utl::getAttribute(child,"type",type,true); + if (type == "navierplate") { + double a = 0.0, b = 0.0, c = 0.0, d = 0.0, t = 0.0; + double E = 10000.0, nu = 0.3, pz = 1.0, xi = 0.0, eta = 0.0; + utl::getAttribute(child,"a",a); + utl::getAttribute(child,"b",b); + utl::getAttribute(child,"c",c); + utl::getAttribute(child,"d",d); + utl::getAttribute(child,"t",t); + utl::getAttribute(child,"E",E); + utl::getAttribute(child,"nu",nu); + utl::getAttribute(child,"pz",pz); + utl::getAttribute(child,"xi",xi); + utl::getAttribute(child,"eta",eta); std::cout <<"\nAnalytic solution: NavierPlate a="<< a <<" b="<< b <<" t="<< t <<" E="<< E <<" nu="<< nu <<" pz="<< pz; - - if (xi != 0.f && eta != 0.f) { + if (xi != 0.0 && eta != 0.0) { std::cout <<" xi="<< xi <<" eta="<< eta; - if (c != 0.f && d != 0.f) { + if (c != 0.0 && d != 0.0) { std::cout <<" c="<< c <<" d="<< d; mySol = new AnaSol(new NavierPlate(a,b,t,E,nu,pz,xi,eta,c,d)); - } else + } + else mySol = new AnaSol(new NavierPlate(a,b,t,E,nu,pz,xi,eta)); - } else + } + else mySol = new AnaSol(new NavierPlate(a,b,t,E,nu,pz)); - } else { - std::cerr <<" ** SIMLinElKL::parse: Unknown analytical solution " - << (child->Attribute("type")?child->Attribute("type"):"") << std::endl; } + else if (type == "expression") { + std::cout <<"\nAnalytical solution: Expression"<< std::endl; + mySol = new AnaSol(child); + } + else + std::cerr <<" ** SIMLinElKL::parse: Unknown analytical solution " + << type <<" (ignored)"<< std::endl; } - else - SIMLinEl2D::parse(child); child = child->NextSiblingElement(); } diff --git a/src/Integrands/AnaSol.C b/src/Integrands/AnaSol.C new file mode 100644 index 00000000..69a61c04 --- /dev/null +++ b/src/Integrands/AnaSol.C @@ -0,0 +1,108 @@ +// $Id$ +//============================================================================== +//! +//! \file AnaSol.C +//! +//! \date Feb 20 2012 +//! +//! \author Knut Morten Okstad / SINTEF +//! +//! \brief Analytical solution fields (primary and secondary). +//! +//============================================================================== + +#include "AnaSol.h" +#include "Functions.h" +#include "Utilities.h" +#include "tinyxml.h" + + +AnaSol::AnaSol (std::istream& is, const int nlines, bool scalarSol) + : scalSol(0), scalSecSol(0), vecSol(0), vecSecSol(0), stressSol(0) +{ + size_t pos = 0; + std::string variables; + for (int i = 0; i < nlines; i++) + { + std::string function = utl::readLine(is); + + if ((pos = function.find("Variables=")) != std::string::npos) + { + variables += function.substr(pos+10); + if (variables[variables.size()-1] != ';') variables += ";"; + std::cout <<"\tVariables="<< variables << std::endl; + } + + if ((pos = function.find("Primary=")) != std::string::npos) + { + std::string primary = function.substr(pos+8); + std::cout <<"\tPrimary="<< primary << std::endl; + if (scalarSol) + scalSol = new EvalFunction((variables+primary).c_str()); + else + vecSol = new VecFuncExpr(primary,variables); + } + + if ((pos = function.find("Secondary=")) != std::string::npos) + { + std::string secondary = function.substr(pos+10); + std::cout <<"\tSecondary="<< secondary << std::endl; + if (scalarSol) + scalSecSol = new VecFuncExpr(secondary,variables); + else + vecSecSol = new TensorFuncExpr(secondary,variables); + } + + if ((pos = function.find("Stress=")) != std::string::npos) + { + std::string stress = function.substr(pos+7); + std::cout <<"\tStress="<< stress << std::endl; + stressSol = new STensorFuncExpr(stress,variables); + } + } +} + + +AnaSol::AnaSol (const TiXmlElement* elem, bool scalarSol) + : scalSol(0), scalSecSol(0), vecSol(0), vecSecSol(0), stressSol(0) +{ + std::string variables; + + const TiXmlElement* var = elem->FirstChildElement("variables"); + if (var && var->FirstChild()) + { + variables = var->FirstChild()->Value(); + if (variables[variables.size()-1] != ';') variables += ";"; + std::cout <<"\tVariables="<< variables << std::endl; + } + + const TiXmlElement* prim = elem->FirstChildElement("primary"); + if (prim && prim->FirstChild()) + { + std::string primary = prim->FirstChild()->Value(); + std::cout <<"\tPrimary="<< primary << std::endl; + if (scalarSol) + scalSol = new EvalFunction((variables+primary).c_str()); + else + vecSol = new VecFuncExpr(primary,variables); + } + + const TiXmlElement* sec = elem->FirstChildElement("secondary"); + if (sec && sec->FirstChild()) + { + std::string secondary = sec->FirstChild()->Value(); + std::cout <<"\tSecondary="<< secondary << std::endl; + if (scalarSol) + scalSecSol = new VecFuncExpr(secondary,variables); + else + vecSecSol = new TensorFuncExpr(secondary,variables); + } + + const TiXmlElement* stress = elem->FirstChildElement("stress"); + if (stress && stress->FirstChild()) + { + std::string stress = sec->FirstChild()->Value(); + std::cout <<"\tStress="<< stress << std::endl; + stressSol = new STensorFuncExpr(stress,variables); + } +} diff --git a/src/Integrands/AnaSol.h b/src/Integrands/AnaSol.h index 9fb0dabc..68e01ba3 100644 --- a/src/Integrands/AnaSol.h +++ b/src/Integrands/AnaSol.h @@ -15,6 +15,9 @@ #define _ANA_SOL_H #include "Function.h" +#include + +class TiXmlElement; /*! @@ -62,6 +65,16 @@ public: AnaSol(STensorFunc* sigma) : scalSol(0), scalSecSol(0), vecSol(0), vecSecSol(0), stressSol(sigma) {} + //! \brief Constructor initializing expression functions by parsing a stream. + //! \param is The file stream to read function definition from + //! \param[in] nlines Number of lines to read + //! \param[in] scalarSol If \e true, the primary solution is a scalar field + AnaSol(std::istream& is, const int nlines, bool scalarSol = true); + //! \brief Constructor initializing expression functions by parsing XML tags. + //! \param[in] elem Pointer to XML-element to extract data from + //! \param[in] scalarSol If \e true, the primary solution is a scalar field + AnaSol(const TiXmlElement* elem, bool scalarSol = true); + //! \brief The destructor frees the analytical solution fields. virtual ~AnaSol() {