diff --git a/dune/stuff/functions.hh b/dune/stuff/functions.hh
index 860bf1d59f60b351a39d20a99ad40ab2c24841e9..73c01129c8ff987a9be9aa55a42430f3c4a54639 100644
--- a/dune/stuff/functions.hh
+++ b/dune/stuff/functions.hh
@@ -17,9 +17,9 @@
#include "functions/checkerboard.hh"
#include "functions/constant.hh"
#include "functions/expression.hh"
-#include "functions/flattop.hh"
#include "functions/ESV2007.hh"
-#include "playground/functions/indicator.hh"
+#include "functions/flattop.hh"
+#include "functions/indicator.hh"
#include "functions/spe10.hh"
namespace Dune {
diff --git a/dune/stuff/functions/ESV2007.hh b/dune/stuff/functions/ESV2007.hh
index e08ef250e4c085860636d60932c63bdb9e3ab40f..2c6a33c1eae6d484707262ae066160f1e6e4f297 100644
--- a/dune/stuff/functions/ESV2007.hh
+++ b/dune/stuff/functions/ESV2007.hh
@@ -10,6 +10,10 @@
#include <type_traits>
+#if HAVE_EIGEN
+#include <Eigen/Eigenvalues>
+#endif
+
#include <dune/geometry/referenceelements.hh>
#include <dune/stuff/common/configuration.hh>
@@ -387,6 +391,230 @@ private:
}; // class Cutoff
+template <class DiffusionFactorType, class DiffusionTensorType>
+class Cutoff
+ : public LocalizableFunctionInterface<typename DiffusionFactorType::EntityType,
+ typename DiffusionFactorType::DomainFieldType, DiffusionFactorType::dimDomain,
+ typename DiffusionFactorType::RangeFieldType, 1, 1>
+{
+ static_assert(std::is_base_of<Tags::LocalizableFunction, DiffusionFactorType>::value,
+ "DiffusionFactorType has to be tagged as a LocalizableFunction!");
+ static_assert(std::is_base_of<Tags::LocalizableFunction, DiffusionTensorType>::value,
+ "DiffusionTensorType has to be tagged as a LocalizableFunction!");
+ typedef typename DiffusionFactorType::EntityType E_;
+ typedef typename DiffusionFactorType::DomainFieldType D_;
+ static const size_t d_ = DiffusionFactorType::dimDomain;
+ typedef typename DiffusionFactorType::RangeFieldType R_;
+ typedef LocalizableFunctionInterface<E_, D_, d_, R_, 1> BaseType;
+ typedef Cutoff<DiffusionFactorType, DiffusionTensorType> ThisType;
+
+ static_assert(DiffusionFactorType::dimRange == 1, "The diffusion factor has to be scalar!");
+ static_assert(DiffusionFactorType::dimRangeCols == 1, "The diffusion factor has to be scalar!");
+
+ static_assert(std::is_same<typename DiffusionTensorType::EntityType, E_>::value, "Types do not match!");
+ static_assert(std::is_same<typename DiffusionTensorType::DomainFieldType, D_>::value, "Types do not match!");
+ static_assert(DiffusionTensorType::dimDomain == d_, "Dimensions do not match!");
+ static_assert(std::is_same<typename DiffusionTensorType::RangeFieldType, R_>::value, "Types do not match!");
+
+ static_assert(DiffusionTensorType::dimRange == d_, "The diffusion tensor has to be a matrix!");
+ static_assert(DiffusionTensorType::dimRangeCols == d_, "The diffusion tensor has to be a matrix!");
+
+ class Localfunction : public LocalfunctionInterface<E_, D_, d_, R_, 1, 1>
+ {
+ typedef LocalfunctionInterface<E_, D_, d_, R_, 1, 1> BaseType;
+
+ public:
+ typedef typename BaseType::EntityType EntityType;
+
+ typedef typename BaseType::DomainFieldType DomainFieldType;
+ static const size_t dimDomain = BaseType::dimDomain;
+ typedef typename BaseType::DomainType DomainType;
+
+ typedef typename BaseType::RangeFieldType RangeFieldType;
+ static const size_t dimRange = BaseType::dimRange;
+ static const size_t dimRangeCols = BaseType::dimRangeCols;
+ typedef typename BaseType::RangeType RangeType;
+
+ typedef typename BaseType::JacobianRangeType JacobianRangeType;
+
+ private:
+ template <class DF, size_t r, size_t rR>
+ struct ComputeDiffusionFactor
+ {
+ static_assert(AlwaysFalse<DF>::value, "Not implemented for these dimensions!");
+ };
+
+ template <class DF>
+ struct ComputeDiffusionFactor<DF, 1, 1>
+ {
+ /**
+ * We try to find the minimum of a polynomial of given order by evaluating it at the points of a quadrature that
+ * would integrate this polynomial exactly.
+ * \todo These are just some heuristics and should be replaced by something proper.
+ */
+ static RangeFieldType min_of(const DF& diffusion_factor, const EntityType& ent)
+ {
+ typename DF::RangeType tmp_value(0);
+ RangeFieldType minimum = std::numeric_limits<RangeFieldType>::max();
+ const auto local_diffusion_factor = diffusion_factor.local_function(ent);
+ const size_t ord = local_diffusion_factor->order();
+ const auto& quadrature =
+ QuadratureRules<DomainFieldType, dimDomain>::rule(ent.type(), boost::numeric_cast<int>(ord));
+ const auto quad_point_it_end = quadrature.end();
+ for (auto quad_point_it = quadrature.begin(); quad_point_it != quad_point_it_end; ++quad_point_it) {
+ local_diffusion_factor->evaluate(quad_point_it->position(), tmp_value);
+ minimum = std::min(minimum, tmp_value[0]);
+ }
+ return minimum;
+ } // ... min_of(...)
+ }; // class ComputeDiffusionFactor< ..., 1, 1 >
+
+ template <class DT, size_t r, size_t rR>
+ struct ComputeDiffusionTensor
+ {
+ static_assert(AlwaysFalse<DT>::value, "Not implemented for these dimensions!");
+ };
+
+ template <class DT, size_t d>
+ struct ComputeDiffusionTensor<DT, d, d>
+ {
+ static RangeFieldType min_eigenvalue_of(const DT& diffusion_tensor, const EntityType& ent)
+ {
+#if !HAVE_EIGEN
+ static_assert(AlwaysFalse<DT>::value, "You are missing eigen!");
+#else
+ const auto local_diffusion_tensor = diffusion_tensor.local_function(ent);
+ assert(local_diffusion_tensor->order() == 0);
+ const auto& reference_element = ReferenceElements<DomainFieldType, dimDomain>::general(ent.type());
+ const Stuff::LA::EigenDenseMatrix<RangeFieldType> tensor =
+ local_diffusion_tensor->evaluate(reference_element.position(0, 0));
+ ::Eigen::EigenSolver<typename Stuff::LA::EigenDenseMatrix<RangeFieldType>::BackendType> eigen_solver(
+ tensor.backend());
+ assert(eigen_solver.info() == ::Eigen::Success);
+ const auto eigenvalues = eigen_solver.eigenvalues(); // <- this should be an Eigen vector of std::complex
+ RangeFieldType min_ev = std::numeric_limits<RangeFieldType>::max();
+ for (size_t ii = 0; ii < boost::numeric_cast<size_t>(eigenvalues.size()); ++ii) {
+ // assert this is real
+ assert(std::abs(eigenvalues[ii].imag()) < 1e-15);
+ // assert that this eigenvalue is positive
+ const RangeFieldType eigenvalue = eigenvalues[ii].real();
+ assert(eigenvalue > 1e-15);
+ min_ev = std::min(min_ev, eigenvalue);
+ }
+ return min_ev;
+#endif // HAVE_EIGEN
+ } // ... min_eigenvalue_of_(...)
+ }; // class Compute< ..., d, d >
+
+ public:
+ Localfunction(const EntityType& ent, const DiffusionFactorType& diffusion_factor,
+ const DiffusionTensorType& diffusion_tensor, const RangeFieldType poincare_constant)
+ : BaseType(ent)
+ , value_(0)
+ {
+ const RangeFieldType min_diffusion_factor =
+ ComputeDiffusionFactor<DiffusionFactorType,
+ DiffusionFactorType::dimRange,
+ DiffusionFactorType::dimRangeCols>::min_of(diffusion_factor, ent);
+ const RangeFieldType min_eigen_value_diffusion_tensor =
+ ComputeDiffusionTensor<DiffusionTensorType,
+ DiffusionTensorType::dimRange,
+ DiffusionTensorType::dimRangeCols>::min_eigenvalue_of(diffusion_tensor, ent);
+ assert(min_diffusion_factor > RangeFieldType(0));
+ assert(min_eigen_value_diffusion_tensor > RangeFieldType(0));
+ const DomainFieldType hh = compute_diameter_of_(ent);
+ value_ = (poincare_constant * hh * hh) / (min_diffusion_factor * min_eigen_value_diffusion_tensor);
+ } // Localfunction(...)
+
+ Localfunction(const Localfunction& /*other*/) = delete;
+
+ Localfunction& operator=(const Localfunction& /*other*/) = delete;
+
+ virtual size_t order() const override final
+ {
+ return 0;
+ }
+
+ virtual void evaluate(const DomainType& UNUSED_UNLESS_DEBUG(xx), RangeType& ret) const override final
+ {
+ assert(this->is_a_valid_point(xx));
+ ret[0] = value_;
+ }
+
+ virtual void jacobian(const DomainType& UNUSED_UNLESS_DEBUG(xx), JacobianRangeType& ret) const override final
+ {
+ assert(this->is_a_valid_point(xx));
+ ret *= RangeFieldType(0);
+ }
+
+ private:
+ static DomainFieldType compute_diameter_of_(const EntityType& ent)
+ {
+ DomainFieldType ret(0);
+ for (auto cc : DSC::valueRange(ent.template count<dimDomain>())) {
+ const auto vertex = ent.template subEntity<dimDomain>(cc)->geometry().center();
+ for (auto dd : DSC::valueRange(cc + 1, ent.template count<dimDomain>())) {
+ const auto other_vertex = ent.template subEntity<dimDomain>(dd)->geometry().center();
+ const auto diff = vertex - other_vertex;
+ ret = std::max(ret, diff.two_norm());
+ }
+ }
+ return ret;
+ } // ... compute_diameter_of_(...)
+
+ RangeFieldType value_;
+ }; // class Localfunction
+
+public:
+ typedef typename BaseType::EntityType EntityType;
+ typedef typename BaseType::LocalfunctionType LocalfunctionType;
+
+ typedef typename BaseType::DomainFieldType DomainFieldType;
+ static const size_t dimDomain = BaseType::dimDomain;
+ typedef typename BaseType::DomainType DomainType;
+
+ typedef typename BaseType::RangeFieldType RangeFieldType;
+ static const size_t dimRange = BaseType::dimRange;
+ static const size_t dimRangeCols = BaseType::dimRangeCols;
+ typedef typename BaseType::RangeType RangeType;
+
+ static std::string static_id()
+ {
+ return BaseType::static_id() + ".ESV2007.cutoff";
+ }
+
+ Cutoff(const DiffusionFactorType& diffusion_factor, const DiffusionTensorType& diffusion_tensor,
+ const RangeFieldType poincare_constant = 1.0 / (M_PIl * M_PIl), const std::string nm = static_id())
+ : diffusion_factor_(diffusion_factor)
+ , diffusion_tensor_(diffusion_tensor)
+ , poincare_constant_(poincare_constant)
+ , name_(nm)
+ {
+ }
+
+ Cutoff(const ThisType& other) = default;
+
+ ThisType& operator=(const ThisType& other) = delete;
+
+ virtual std::string name() const override final
+ {
+ return name_;
+ }
+
+ virtual std::unique_ptr<LocalfunctionType> local_function(const EntityType& entity) const override final
+ {
+ return std::unique_ptr<Localfunction>(
+ new Localfunction(entity, diffusion_factor_, diffusion_tensor_, poincare_constant_));
+ }
+
+private:
+ const DiffusionFactorType& diffusion_factor_;
+ const DiffusionTensorType& diffusion_tensor_;
+ const RangeFieldType poincare_constant_;
+ std::string name_;
+}; // class Cutoff
+
+
} // namespace ESV2007
} // namespace Functions
} // namespace Stuff
diff --git a/dune/stuff/functions/expression.hh b/dune/stuff/functions/expression.hh
index ed67d39f67dc9d93530b9dc339ba120e621c98fc..ebfd6c2b15fb81de8f9c115771b267da9cbaee1c 100644
--- a/dune/stuff/functions/expression.hh
+++ b/dune/stuff/functions/expression.hh
@@ -214,24 +214,24 @@ public:
#ifndef NDEBUG
#ifndef DUNE_STUFF_FUNCTIONS_EXPRESSION_DISABLE_CHECKS
bool failure = false;
- std::string type;
+ std::string error_type;
for (size_t rr = 0; rr < dimRange; ++rr) {
tmp_row_ = ret[rr];
for (size_t cc = 0; cc < dimRangeCols; ++cc) {
if (DSC::isnan(tmp_row_[cc])) {
- failure = true;
- type = "NaN";
+ failure = true;
+ error_type = "NaN";
} else if (DSC::isinf(tmp_row_[cc])) {
- failure = true;
- type = "inf";
+ failure = true;
+ error_type = "inf";
} else if (std::abs(tmp_row_[cc]) > (0.9 * std::numeric_limits<double>::max())) {
- failure = true;
- type = "an unlikely value";
+ failure = true;
+ error_type = "an unlikely value";
}
if (failure)
DUNE_THROW(Stuff::Exceptions::internal_error,
"evaluating this function yielded "
- << type
+ << error_type
<< "!\n"
<< "The variable of this function is: "
<< function_->variable()
diff --git a/dune/stuff/functions/expression/mathexpr.cc b/dune/stuff/functions/expression/mathexpr.cc
index a099475e187ca5e89c3d0e76c2a7f384f9fe8c42..79fbede1b9badb16bf2fed48e147206e69dcab1c 100644
--- a/dune/stuff/functions/expression/mathexpr.cc
+++ b/dune/stuff/functions/expression/mathexpr.cc
@@ -21,622 +21,336 @@ This software comes with absolutely no warranty.
#include "mathexpr.hh"
-char* MidStr(const char* s, int i1, int i2)
+// clang-format off
+
+char* MidStr(const char*s,int i1,int i2)
{
- if (i1 < 0 || i2 >= (int)strlen(s) || i1 > i2) {
- char* cp = new char[1];
- cp[0] = '\0';
- return cp;
- }
- char* s1 = new char[i2 - i1 + 2];
- int i;
- for (i = i1; i <= i2; i++)
- s1[i - i1] = s[i];
- s1[i2 - i1 + 1] = 0;
- return s1;
+ if(i1<0||i2>=(int)strlen(s)||i1>i2){
+ char* cp = new char[1];
+ cp[0] = '\0';
+ return cp;
+ }
+ char*s1=new char[i2-i1+2];
+ int i;
+ for(i=i1;i<=i2;i++)s1[i-i1]=s[i];
+ s1[i2-i1+1]=0;return s1;
}
-char* CopyStr(const char* s)
-{
- char* s1 = new char[strlen(s) + 1];
- char* s12 = s1;
- const char* s2 = s;
- while ((*s12++ = *s2++))
+char* CopyStr(const char*s)
+{char*s1=new char[strlen(s)+1];char*s12=s1;const char*s2=s;
+while((*s12++=*s2++))
;
- return s1;
-}
+return s1;}
-void InsStr(char*& s, int n, char c) // Warning : deletes the old string
-{
- if (n < 0 || n > (int)strlen(s))
- return;
- char* s1 = new char[strlen(s) + 2];
- int i;
- for (i = 0; i < n; i++)
- s1[i] = s[i];
- s1[n] = c;
- for (i = n + 1; s[i - 1]; i++)
- s1[i] = s[i - 1];
- s1[i] = 0;
- delete[] s;
- s = s1;
+void InsStr(char*&s,int n,char c)// Warning : deletes the old string
+{if(n<0||n>(int)strlen(s))return;
+char*s1=new char[strlen(s)+2];
+int i;
+for(i=0;i<n;i++)s1[i]=s[i];
+s1[n]=c;for(i=n+1;s[i-1];i++)s1[i]=s[i-1];
+s1[i]=0;
+delete[]s;s=s1;
}
-signed char EqStr(const char* s, const char* s2)
-{
- if (strlen(s) != strlen(s2))
- return 0;
- int i;
- for (i = 0; s[i]; i++)
- if (s[i] != s2[i])
- return 0;
- return 1;
+signed char EqStr(const char*s,const char*s2)
+{if(strlen(s)!=strlen(s2))return 0;
+int i;
+for(i=0;s[i];i++)if(s[i]!=s2[i])return 0;
+return 1;
}
-signed char CompStr(const char* s, int n, const char* s2)
-{
- if (n < 0 || n >= (int)strlen(s) || n + (int)strlen(s2) > (int)strlen(s))
- return 0;
- int i;
- for (i = 0; s2[i]; i++)
- if (s[i + n] != s2[i])
- return 0;
- return 1;
+signed char CompStr(const char*s,int n,const char*s2)
+{if(n<0||n>=(int)strlen(s)||n+(int)strlen(s2)>(int)strlen(s))return 0;
+int i;
+for(i=0;s2[i];i++)if(s[i+n]!=s2[i])return 0;
+return 1;
}
-void DelStr(char*& s, int n) // Deletes the old string
-{
- char* s1 = new char[strlen(s)];
- int i;
- for (i = 0; i < n; i++)
- s1[i] = s[i];
- for (i = n; s[i + 1]; i++)
- s1[i] = s[i + 1];
- s1[i] = 0;
- delete[] s;
- s = s1;
+void DelStr(char*&s,int n)//Deletes the old string
+{char*s1=new char[strlen(s)];
+int i;
+for(i=0;i<n;i++)s1[i]=s[i];
+for(i=n;s[i+1];i++)s1[i]=s[i+1];
+s1[i]=0;
+delete[]s;s=s1;
}
-RVar::RVar(const RVar& rvarp)
-{
- if (this == &rvarp)
- return;
- pval = rvarp.pval;
- name = CopyStr(rvarp.name);
+RVar::RVar(const RVar & rvarp)
+{if(this==&rvarp)return;pval=rvarp.pval;name=CopyStr(rvarp.name);
}
-RVar::RVar(const char* namep, double* pvalp)
-{
- pval = pvalp;
- name = CopyStr(namep);
-}
+RVar::RVar(const char*namep,double*pvalp)
+{pval=pvalp;name=CopyStr(namep);}
RVar::~RVar()
-{
- if (name != NULL)
- delete[] name;
-}
+{if(name!=NULL)delete[] name;}
RFunction::RFunction()
{
- type = -1;
- name = new char[1];
- name[0] = 0;
- nvars = 0;
- ppvar = NULL;
- pfuncval = NULL;
- op = ErrVal;
- buf = NULL;
+ type=-1;name=new char[1];name[0]=0;
+ nvars=0;ppvar=NULL;pfuncval=NULL;op=ErrVal;buf=NULL;
}
-RFunction::RFunction(double((*pfuncvalp)(double)))
+RFunction::RFunction(double ((*pfuncvalp)(double)))
{
- type = 0;
- pfuncval = pfuncvalp;
- name = new char[1];
- name[0] = 0;
- nvars = 1;
- ppvar = NULL;
- op = ErrVal;
- buf = NULL;
+ type=0;pfuncval=pfuncvalp;name=new char[1];name[0]=0;
+ nvars=1;ppvar=NULL;op=ErrVal;buf=NULL;
}
RFunction::RFunction(const RFunction& rfunc)
{
- if (this == &rfunc)
- return;
- type = rfunc.type;
- op = rfunc.op;
- pfuncval = rfunc.pfuncval;
- name = CopyStr(rfunc.name);
- nvars = rfunc.nvars;
- if (rfunc.ppvar != NULL && nvars) {
- ppvar = new PRVar[nvars];
- int i;
- for (i = 0; i < nvars; i++)
- ppvar[i] = rfunc.ppvar[i];
- buf = new double[nvars];
- } else {
- ppvar = NULL;
- buf = NULL;
- }
+ if(this==&rfunc)return;
+ type=rfunc.type;op=rfunc.op;
+ pfuncval=rfunc.pfuncval;
+ name=CopyStr(rfunc.name);
+ nvars=rfunc.nvars;
+ if(rfunc.ppvar!=NULL&&nvars){
+ ppvar=new PRVar[nvars];
+ int i;for(i=0;i<nvars;i++)ppvar[i]=rfunc.ppvar[i];
+ buf=new double[nvars];
+ }else {ppvar=NULL;buf=NULL;}
}
-RFunction::RFunction(const ROperation& opp, RVar* pvarp)
- : op(opp)
+RFunction::RFunction(const ROperation& opp,RVar* pvarp):op(opp)
{
- type = 1;
- name = new char[1];
- name[0] = 0;
- nvars = 1;
- ppvar = new PRVar[1];
- ppvar[0] = pvarp;
- buf = new double[1];
+ type=1;name=new char[1];name[0]=0;
+ nvars=1;ppvar=new PRVar[1];ppvar[0]=pvarp;buf=new double[1];
}
-RFunction::RFunction(const ROperation& opp, int nvarsp, RVar** ppvarp)
- : op(opp)
+RFunction::RFunction(const ROperation& opp, int nvarsp,RVar**ppvarp):op(opp)
{
- type = 1;
- name = new char[1];
- name[0] = 0;
- nvars = nvarsp;
- if (nvars) {
- ppvar = new PRVar[nvars];
- int i;
- for (i = 0; i < nvars; i++)
- ppvar[i] = ppvarp[i];
- buf = new double[nvars];
- } else {
- ppvar = NULL;
- buf = NULL;
- }
+ type=1;name=new char[1];name[0]=0;
+ nvars=nvarsp;
+ if(nvars){
+ ppvar=new PRVar[nvars];
+ int i;for(i=0;i<nvars;i++)ppvar[i]=ppvarp[i];
+ buf=new double[nvars];
+ }else {ppvar=NULL;buf=NULL;}
}
RFunction::~RFunction()
{
- if (name != NULL)
- delete[] name;
- if (ppvar != NULL)
- delete[] ppvar;
- if (buf != NULL)
- delete[] buf;
+ if(name!=NULL)delete[]name;
+ if(ppvar!=NULL)delete[]ppvar;
+ if(buf!=NULL)delete[]buf;
}
RFunction& RFunction::operator=(const RFunction& rfunc)
{
- if (this == &rfunc)
- return *this;
- type = rfunc.type;
- op = rfunc.op;
- pfuncval = rfunc.pfuncval;
- delete[] name;
- name = CopyStr(rfunc.name);
- if (ppvar != NULL)
- delete[] ppvar;
- ppvar = NULL;
- if (buf != NULL)
- delete[] buf;
- buf = NULL;
- nvars = rfunc.nvars;
- if (type == 1 && nvars) {
- ppvar = new PRVar[nvars];
- buf = new double[nvars];
- int i;
- for (i = 0; i < nvars; i++)
- ppvar[i] = rfunc.ppvar[i];
+ if(this==&rfunc)return *this;
+ type=rfunc.type;op=rfunc.op;
+ pfuncval=rfunc.pfuncval;
+ delete[]name;
+ name=CopyStr(rfunc.name);
+ if(ppvar!=NULL)delete[]ppvar;
+ ppvar=NULL;
+ if(buf!=NULL)delete[]buf;
+ buf=NULL;
+ nvars=rfunc.nvars;
+ if(type==1&&nvars){
+ ppvar=new PRVar[nvars];buf=new double[nvars];
+ int i;for(i=0;i<nvars;i++)ppvar[i]=rfunc.ppvar[i];
}
return *this;
}
-void RFunction::SetName(const char* s)
-{
- if (name != NULL)
- delete[] name;
- name = CopyStr(s);
-}
+void RFunction::SetName(const char*s)
+{if(name!=NULL)delete[]name;name=CopyStr(s);}
double RFunction::Val(double x) const
{
- if (type == -1 || nvars >= 2)
- return ErrVal;
- if (type == 0)
- return (*pfuncval)(x);
- double xb = *(*ppvar)->pval, y;
- *(*ppvar)->pval = x; // Warning : could cause trouble if this value is used in a parallel process
- y = op.Val();
- *(*ppvar)->pval = xb;
+ if(type==-1||nvars>=2)return ErrVal;
+ if(type==0)return (*pfuncval)(x);
+ double xb=*(*ppvar)->pval,y;
+ *(*ppvar)->pval=x; // Warning : could cause trouble if this value is used in a parallel process
+ y=op.Val();
+ *(*ppvar)->pval=xb;
return y;
}
-double RFunction::Val(double* pv) const
+double RFunction::Val(double*pv) const
{
- if (type == -1)
- return ErrVal;
- if (type == 0)
- return (*pfuncval)(*pv);
+ if(type==-1)return ErrVal;
+ if(type==0)return (*pfuncval)(*pv);
double y;
int i;
- for (i = 0; i < nvars; i++) {
- buf[i] = *ppvar[i]->pval;
+ for(i=0;i<nvars;i++){
+ buf[i]=*ppvar[i]->pval;
// Warning : could cause trouble if this value is used in a parallel process
- *ppvar[i]->pval = pv[i];
+ *ppvar[i]->pval=pv[i];
}
- y = op.Val();
- for (i = 0; i < nvars; i++)
- *ppvar[i]->pval = buf[i];
+ y=op.Val();
+ for(i=0;i<nvars;i++)*ppvar[i]->pval=buf[i];
return y;
}
ROperation::ROperation()
-{
- op = ErrOp;
- mmb1 = NULL;
- mmb2 = NULL;
- ValC = ErrVal;
- pvar = NULL;
- pvarval = NULL;
- pfunc = NULL;
- containfuncflag = 0;
- pinstr = NULL;
- pvals = NULL;
- ppile = NULL;
- pfuncpile = NULL;
- BuildCode();
-}
+{op=ErrOp;mmb1=NULL;mmb2=NULL;ValC=ErrVal;pvar=NULL;pvarval=NULL;pfunc=NULL;containfuncflag=0;pinstr=NULL;pvals=NULL;ppile=NULL;pfuncpile=NULL;BuildCode();}
ROperation::~ROperation()
{
Destroy();
}
-ROperation::ROperation(const ROperation& ROp)
+ROperation::ROperation(const ROperation&ROp)
{
- op = ROp.op;
- pvar = ROp.pvar;
- pvarval = ROp.pvarval;
- ValC = ROp.ValC;
- pfunc = ROp.pfunc;
- containfuncflag = 0;
- pinstr = NULL;
- pvals = NULL;
- ppile = NULL;
- pfuncpile = NULL;
- if (ROp.mmb1 != NULL)
- mmb1 = new ROperation(*(ROp.mmb1));
- else
- mmb1 = NULL;
- if (ROp.mmb2 != NULL)
- mmb2 = new ROperation(*(ROp.mmb2));
- else
- mmb2 = NULL;
+ op=ROp.op;pvar=ROp.pvar;pvarval=ROp.pvarval;ValC=ROp.ValC;pfunc=ROp.pfunc;containfuncflag=0;pinstr=NULL;pvals=NULL;ppile=NULL;pfuncpile=NULL;
+ if(ROp.mmb1!=NULL)mmb1=new ROperation(*(ROp.mmb1));else mmb1=NULL;
+ if(ROp.mmb2!=NULL)mmb2=new ROperation(*(ROp.mmb2));else mmb2=NULL;
BuildCode();
}
ROperation::ROperation(double x)
{
- if (x == ErrVal) {
- op = ErrOp;
- mmb1 = NULL;
- mmb2 = NULL;
- ValC = ErrVal;
- } else if (x >= 0) {
- op = Num;
- mmb1 = NULL;
- mmb2 = NULL;
- ValC = x;
- } else {
- op = Opp;
- mmb1 = NULL;
- mmb2 = new ROperation(-x);
- ValC = ErrVal;
- }
- pvar = NULL;
- pvarval = NULL;
- pfunc = NULL;
- containfuncflag = 0;
- pinstr = NULL;
- pvals = NULL;
- ppile = NULL;
- pfuncpile = NULL;
+ if(x==ErrVal){op=ErrOp;mmb1=NULL;mmb2=NULL;ValC=ErrVal;}
+ else if(x>=0){op=Num;mmb1=NULL;mmb2=NULL;ValC=x;}
+ else{op=Opp;mmb1=NULL;mmb2=new ROperation(-x);ValC=ErrVal;}
+ pvar=NULL;pvarval=NULL;pfunc=NULL;containfuncflag=0;pinstr=NULL;pvals=NULL;ppile=NULL;pfuncpile=NULL;
BuildCode();
}
-ROperation::ROperation(const RVar& varp)
-{
- op = Var;
- mmb1 = NULL;
- mmb2 = NULL;
- ValC = ErrVal;
- pvar = &varp;
- pvarval = varp.pval;
- containfuncflag = 0;
- pfunc = NULL;
- pinstr = NULL;
- pvals = NULL;
- ppile = NULL;
- pfuncpile = NULL;
- BuildCode();
-}
+ROperation::ROperation(const RVar&varp)
+{op=Var;mmb1=NULL;mmb2=NULL;ValC=ErrVal;pvar=&varp;pvarval=varp.pval;containfuncflag=0;pfunc=NULL;pinstr=NULL;pvals=NULL;ppile=NULL;pfuncpile=NULL;BuildCode();}
ROperation& ROperation::operator=(const ROperation& ROp)
{
- if (this == &ROp)
- return *this;
+ if(this==&ROp)return *this;
Destroy();
- op = ROp.op;
- pvar = ROp.pvar;
- pvarval = ROp.pvarval;
- ValC = ROp.ValC;
- pfunc = ROp.pfunc;
- containfuncflag = 0;
- pinstr = NULL;
- pvals = NULL;
- ppile = NULL;
- pfuncpile = NULL;
- if (ROp.mmb1 != NULL)
- mmb1 = new ROperation(*(ROp.mmb1));
- else
- mmb1 = NULL;
- if (ROp.mmb2 != NULL)
- mmb2 = new ROperation(*(ROp.mmb2));
- else
- mmb2 = NULL;
+ op=ROp.op;pvar=ROp.pvar;pvarval=ROp.pvarval;ValC=ROp.ValC;pfunc=ROp.pfunc;containfuncflag=0;pinstr=NULL;pvals=NULL;ppile=NULL;pfuncpile=NULL;
+ if(ROp.mmb1!=NULL)mmb1=new ROperation(*(ROp.mmb1));else mmb1=NULL;
+ if(ROp.mmb2!=NULL)mmb2=new ROperation(*(ROp.mmb2));else mmb2=NULL;
BuildCode();
return *this;
}
-int operator==(const ROperation& op, const double v)
-{
- return (op.op == Num && op.ValC == v);
+int operator==(const ROperation& op,const double v)
+{return(op.op==Num&&op.ValC==v);}
+
+int operator==(const ROperation& op1,const ROperation& op2)
+{if(op1.op!=op2.op)return 0;
+if(op1.op==Var)return(*(op1.pvar)==*(op2.pvar));
+if(op1.op==Fun)return(op1.pfunc==op2.pfunc); // *op1.pfunc==*op2.pfunc could imply infinite loops in cases of self-dependence
+if(op1.op==Num)return(op1.ValC==op2.ValC);
+if(op1.mmb1==NULL&&op2.mmb1!=NULL)return 0;
+if(op1.mmb2==NULL&&op2.mmb2!=NULL)return 0;
+if(op2.mmb1==NULL&&op1.mmb1!=NULL)return 0;
+if(op2.mmb2==NULL&&op1.mmb2!=NULL)return 0;
+return(((op1.mmb1==NULL&&op2.mmb1==NULL)||(*(op1.mmb1)==*(op2.mmb1)))&&
+ ((op1.mmb2==NULL&&op2.mmb2==NULL)||(*(op1.mmb2)==*(op2.mmb2))));
}
-int operator==(const ROperation& op1, const ROperation& op2)
-{
- if (op1.op != op2.op)
- return 0;
- if (op1.op == Var)
- return (*(op1.pvar) == *(op2.pvar));
- if (op1.op == Fun)
- return (op1.pfunc == op2.pfunc); // *op1.pfunc==*op2.pfunc could imply infinite loops in cases of self-dependence
- if (op1.op == Num)
- return (op1.ValC == op2.ValC);
- if (op1.mmb1 == NULL && op2.mmb1 != NULL)
- return 0;
- if (op1.mmb2 == NULL && op2.mmb2 != NULL)
- return 0;
- if (op2.mmb1 == NULL && op1.mmb1 != NULL)
- return 0;
- if (op2.mmb2 == NULL && op1.mmb2 != NULL)
- return 0;
- return (((op1.mmb1 == NULL && op2.mmb1 == NULL) || (*(op1.mmb1) == *(op2.mmb1)))
- && ((op1.mmb2 == NULL && op2.mmb2 == NULL) || (*(op1.mmb2) == *(op2.mmb2))));
-}
-
-int operator!=(const ROperation& op1, const ROperation& op2)
+int operator!=(const ROperation& op1,const ROperation& op2)
{
- if (op1.op != op2.op)
- return 1;
- if (op1.op == Var)
- return (op1.pvar != op2.pvar);
- if (op1.op == Fun)
- return (!(op1.pfunc == op2.pfunc)); // *op1.pfunc==*op2.pfunc could imply infinite loops in cases of self-dependence
- if (op1.op == Num)
- return (op1.ValC != op2.ValC);
- if (op1.mmb1 == NULL && op2.mmb1 != NULL)
- return 1;
- if (op1.mmb2 == NULL && op2.mmb2 != NULL)
- return 1;
- if (op2.mmb1 == NULL && op1.mmb1 != NULL)
- return 1;
- if (op2.mmb2 == NULL && op1.mmb2 != NULL)
- return 1;
- return (((op1.mmb1 != NULL || op2.mmb1 != NULL) && (*(op1.mmb1) != *(op2.mmb1)))
- || ((op1.mmb2 != NULL || op2.mmb2 != NULL) && (*(op1.mmb2) != *(op2.mmb2))));
+ if(op1.op!=op2.op)return 1;
+ if(op1.op==Var)return(op1.pvar!=op2.pvar);
+ if(op1.op==Fun)return(!(op1.pfunc==op2.pfunc)); // *op1.pfunc==*op2.pfunc could imply infinite loops in cases of self-dependence
+ if(op1.op==Num)return(op1.ValC!=op2.ValC);
+ if(op1.mmb1==NULL&&op2.mmb1!=NULL)return 1;
+ if(op1.mmb2==NULL&&op2.mmb2!=NULL)return 1;
+ if(op2.mmb1==NULL&&op1.mmb1!=NULL)return 1;
+ if(op2.mmb2==NULL&&op1.mmb2!=NULL)return 1;
+ return(((op1.mmb1!=NULL||op2.mmb1!=NULL)&&(*(op1.mmb1)!=*(op2.mmb1)))||
+ ((op1.mmb2!=NULL||op2.mmb2!=NULL)&&(*(op1.mmb2)!=*(op2.mmb2))));
}
ROperation ROperation::operator+() const
-{
- return *this;
-}
+{return *this;}
ROperation ROperation::operator-() const
-{
- if (op == Num)
- return -ValC;
- ROperation resultat;
- if (op == Opp)
- resultat = *mmb2;
- else {
- resultat.op = Opp;
- resultat.mmb2 = new ROperation(*this);
- };
- return resultat;
+{if(op==Num)return -ValC;
+ROperation resultat;
+if(op==Opp)resultat=*mmb2;else{resultat.op=Opp;resultat.mmb2=new ROperation(*this);};
+return resultat;
}
-ROperation operator, (const ROperation& op1, const ROperation& op2)
-{
- ROperation resultat;
- resultat.op = Juxt;
- resultat.mmb1 = new ROperation(op1);
- resultat.mmb2 = new ROperation(op2);
- return resultat;
+ROperation operator,(const ROperation& op1,const ROperation& op2)
+{ROperation resultat;
+resultat.op=Juxt;resultat.mmb1=new ROperation(op1);
+resultat.mmb2=new ROperation(op2);
+return resultat;
}
-ROperation operator+(const ROperation& op1, const ROperation& op2)
-{
- if (op1.op == Num && op2.op == Num)
- return op1.ValC + op2.ValC;
- if (op1 == 0.)
- return op2;
- if (op2 == 0.)
- return op1;
- if (op1.op == Opp)
- return op2 - *(op1.mmb2);
- if (op2.op == Opp)
- return op1 - *(op2.mmb2);
- ROperation resultat;
- resultat.op = Add;
- resultat.mmb1 = new ROperation(op1);
- resultat.mmb2 = new ROperation(op2);
- return resultat;
-}
-
-ROperation operator-(const ROperation& op1, const ROperation& op2)
-{
- if (op1.op == Num && op2.op == Num)
- return op1.ValC - op2.ValC;
- if (op1 == 0.)
- return -op2;
- if (op2 == 0.)
- return op1;
- if (op1.op == Opp)
- return -(op2 + *(op1.mmb2));
- if (op2.op == Opp)
- return op1 + *(op2.mmb2);
- ROperation resultat;
- resultat.op = Sub;
- resultat.mmb1 = new ROperation(op1);
- resultat.mmb2 = new ROperation(op2);
- return resultat;
-}
-
-ROperation operator*(const ROperation& op1, const ROperation& op2)
-{
- if (op1.op == Num && op2.op == Num)
- return op1.ValC * op2.ValC;
- if (op1 == 0. || op2 == 0.)
- return 0.;
- if (op1 == 1.)
- return op2;
- if (op2 == 1.)
- return op1;
- if (op1.op == Opp)
- return -(*(op1.mmb2) * op2);
- if (op2.op == Opp)
- return -(op1 * *(op2.mmb2));
- ROperation resultat;
- resultat.op = Mult;
- resultat.mmb1 = new ROperation(op1);
- resultat.mmb2 = new ROperation(op2);
- return resultat;
-}
-
-ROperation operator/(const ROperation& op1, const ROperation& op2)
+ROperation operator+(const ROperation& op1,const ROperation& op2)
{
- if (op1.op == Num && op2.op == Num)
- return (op2.ValC ? op1.ValC / op2.ValC : ErrVal);
- if (op1 == 0.0)
- return 0.;
- if (op2 == 1.)
- return op1;
- if (op2 == 0.)
- return ErrVal;
- if (op1.op == Opp)
- return -(*(op1.mmb2) / op2);
- if (op2.op == Opp)
- return -(op1 / (*(op2.mmb2)));
- ROperation resultat;
- resultat.op = Div;
- resultat.mmb1 = new ROperation(op1);
- resultat.mmb2 = new ROperation(op2);
- return resultat;
-}
-
-ROperation operator^(const ROperation& op1, const ROperation& op2) {
- if (op1 == 0.)
- return 0.;
- if (op2 == 0.)
- return 1.;
- if (op2 == 1.)
- return op1;
- ROperation resultat;
- resultat.op = Pow;
- resultat.mmb1 = new ROperation(op1);
- resultat.mmb2 = new ROperation(op2);
- return resultat;
+if(op1.op==Num&&op2.op==Num)return op1.ValC+op2.ValC;
+if(op1==0.)return op2;if(op2==0.)return op1;
+if(op1.op==Opp)return op2-*(op1.mmb2);if(op2.op==Opp)return op1-*(op2.mmb2);
+ROperation resultat;
+resultat.op=Add;resultat.mmb1=new ROperation(op1);
+resultat.mmb2=new ROperation(op2);
+return resultat;
}
-ROperation sqrt(const ROperation& op)
+ROperation operator-(const ROperation& op1,const ROperation& op2)
{
- ROperation rop;
- rop.op = Sqrt;
- rop.mmb2 = new ROperation(op);
- return rop;
+if(op1.op==Num&&op2.op==Num)return op1.ValC-op2.ValC;
+if(op1==0.)return -op2;if(op2==0.)return op1;
+if(op1.op==Opp)return -(op2+*(op1.mmb2));if(op2.op==Opp)return op1+*(op2.mmb2);
+ROperation resultat;
+resultat.op=Sub;resultat.mmb1=new ROperation(op1);
+resultat.mmb2=new ROperation(op2);
+return resultat;
}
-ROperation abs(const ROperation& op)
+
+ROperation operator*(const ROperation& op1,const ROperation& op2)
{
- ROperation rop;
- rop.op = Abs;
- rop.mmb2 = new ROperation(op);
- return rop;
+if(op1.op==Num&&op2.op==Num)return op1.ValC*op2.ValC;
+if(op1==0.||op2==0.)return 0.;
+if(op1==1.)return op2;if(op2==1.)return op1;
+if(op1.op==Opp)return -(*(op1.mmb2)*op2);if(op2.op==Opp)return -(op1**(op2.mmb2));
+ROperation resultat;
+resultat.op=Mult;resultat.mmb1=new ROperation(op1);
+resultat.mmb2=new ROperation(op2);
+return resultat;
}
-ROperation sin(const ROperation& op)
-{
- ROperation rop;
- rop.op = Sin;
- rop.mmb2 = new ROperation(op);
- return rop;
+
+ROperation operator/(const ROperation& op1,const ROperation& op2)
+{if(op1.op==Num&&op2.op==Num)return (op2.ValC?op1.ValC/op2.ValC:ErrVal);
+if(op1==0.0)return 0.;if(op2==1.)return op1;if(op2==0.)return ErrVal;
+if(op1.op==Opp)return -(*(op1.mmb2)/op2);if(op2.op==Opp)return -(op1/(*(op2.mmb2)));
+ROperation resultat;
+resultat.op=Div;resultat.mmb1=new ROperation(op1);
+resultat.mmb2=new ROperation(op2);
+return resultat;
}
-ROperation cos(const ROperation& op)
-{
- ROperation rop;
- rop.op = Cos;
- rop.mmb2 = new ROperation(op);
- return rop;
+
+ROperation operator^(const ROperation& op1,const ROperation& op2)
+{if(op1==0.)return 0.;
+if(op2==0.)return 1.;
+if(op2==1.)return op1;
+ROperation resultat;
+resultat.op=Pow;resultat.mmb1=new ROperation(op1);
+resultat.mmb2=new ROperation(op2);
+return resultat;
}
+
+ROperation sqrt(const ROperation& op)
+{ROperation rop;rop.op=Sqrt;rop.mmb2=new ROperation(op);return rop;}
+ROperation abs(const ROperation& op)
+{ROperation rop;rop.op=Abs;rop.mmb2=new ROperation(op);return rop;}
+ROperation sin(const ROperation& op)
+{ROperation rop;rop.op=Sin;rop.mmb2=new ROperation(op);return rop;}
+ROperation cos(const ROperation& op)
+{ROperation rop;rop.op=Cos;rop.mmb2=new ROperation(op);return rop;}
ROperation tan(const ROperation& op)
-{
- ROperation rop;
- rop.op = Tg;
- rop.mmb2 = new ROperation(op);
- return rop;
-}
+{ROperation rop;rop.op=Tg;rop.mmb2=new ROperation(op);return rop;}
ROperation log(const ROperation& op)
-{
- ROperation rop;
- rop.op = Ln;
- rop.mmb2 = new ROperation(op);
- return rop;
-}
+{ROperation rop;rop.op=Ln;rop.mmb2=new ROperation(op);return rop;}
ROperation exp(const ROperation& op)
-{
- ROperation rop;
- rop.op = Exp;
- rop.mmb2 = new ROperation(op);
- return rop;
-}
+{ROperation rop;rop.op=Exp;rop.mmb2=new ROperation(op);return rop;}
ROperation acos(const ROperation& op)
-{
- ROperation rop;
- rop.op = Acos;
- rop.mmb2 = new ROperation(op);
- return rop;
-}
+{ROperation rop;rop.op=Acos;rop.mmb2=new ROperation(op);return rop;}
ROperation asin(const ROperation& op)
-{
- ROperation rop;
- rop.op = Asin;
- rop.mmb2 = new ROperation(op);
- return rop;
-}
+{ROperation rop;rop.op=Asin;rop.mmb2=new ROperation(op);return rop;}
ROperation atan(const ROperation& op)
-{
- ROperation rop;
- rop.op = Atan;
- rop.mmb2 = new ROperation(op);
- return rop;
-}
+{ROperation rop;rop.op=Atan;rop.mmb2=new ROperation(op);return rop;}
-ROperation ApplyOperator(int n, ROperation** pops, ROperation (*func)(const ROperation&, const ROperation&))
+ROperation ApplyOperator(int n,ROperation**pops,ROperation (*func)(const ROperation&,const ROperation&))
{
- if (n <= 0)
- return ErrVal;
- if (n == 1)
- return *pops[0];
- if (n == 2)
- return (*func)(*pops[0], *pops[1]);
- return (*func)(*pops[0], ApplyOperator(n - 1, pops + 1, func));
+ if(n<=0)return ErrVal;
+ if(n==1)return *pops[0];
+ if(n==2)return (*func)(*pops[0],*pops[1]);
+ return (*func)(*pops[0],ApplyOperator(n-1,pops+1,func));
}
ROperation RFunction::operator()(const ROperation& _op)
@@ -657,251 +371,151 @@ ROperation RFunction::operator()(const ROperation& _op)
delete[]ppvar2;
return op2;
*/
- ROperation op2;
- op2.op = Fun;
- op2.pfunc = this;
- op2.mmb2 = new ROperation(_op);
+ ROperation op2;op2.op=Fun;op2.pfunc=this;op2.mmb2=new ROperation(_op);
return op2;
}
-// Auxiliary string functions
+//Auxiliary string functions
-void SupprSpaces(char*& s) // Deletes the old string
+void SupprSpaces(char*&s)//Deletes the old string
{
- int i;
- for (i = 0; s[i]; i++)
- if (s[i] == ' ' || s[i] == '\t' || s[i] == '\n')
- DelStr(s, i--);
+ int i;
+ for(i=0;s[i];i++)if(s[i]==' '||s[i]=='\t'||s[i]=='\n')DelStr(s,i--);
}
signed char IsNumeric(char c)
-{
- if (c != '0' && c != '1' && c != '2' && c != '3' && c != '4' && c != '5' && c != '6' && c != '7' && c != '8'
- && c != '9'
- && c != '.')
- return 0;
- return 1;
+{if(c!='0'&&c!='1'&&c!='2'&&c!='3'&&c!='4'
+ &&c!='5'&&c!='6'&&c!='7'&&c!='8'&&c!='9'&&c!='.')return 0;
+return 1;
}
-signed char IsTNumeric(char* s)
-{
- int i;
- for (i = 0; i < (int)strlen(s); i++)
- if (!IsNumeric(s[i]))
- return 0;
- return 1;
+signed char IsTNumeric(char *s)
+{int i;for(i=0;i<(int)strlen(s);i++)if(!IsNumeric(s[i]))return 0;return 1;
}
-int SearchCorOpenbracket(char* s, int n) // Searchs the corresponding bracket of an opening bracket
-{
- if (n >= (int)strlen(s) - 1)
- return -1;
- int i, c = 1;
- for (i = n + 1; s[i]; i++) {
- if (s[i] == '(')
- c++;
- else if (s[i] == ')')
- c--;
- if (!c)
- return i;
+int SearchCorOpenbracket(char*s,int n) //Searchs the corresponding bracket of an opening bracket
+{if(n>=(int)strlen(s)-1)return -1;
+int i,c=1;
+for(i=n+1;s[i];i++){
+ if(s[i]=='(')c++;else if(s[i]==')')c--;
+ if(!c)return i;
};
- return -1;
+return -1;
}
-int SearchCorClosebracket(char* s, int n) // Searchs the corresponding bracket of a closing bracket
-{
- if (n < 1)
- return -1;
- int i, c = 1;
- for (i = n - 1; i >= 0; i--) {
- if (s[i] == ')')
- c++;
- else if (s[i] == '(')
- c--;
- if (!c)
- return i;
+int SearchCorClosebracket(char*s,int n) //Searchs the corresponding bracket of a closing bracket
+{if(n<1)return -1;
+int i,c=1;
+for(i=n-1;i>=0;i--){
+ if(s[i]==')')c++;else if(s[i]=='(')c--;
+ if(!c)return i;
};
- return -1;
+return -1;
}
-int SearchOperator(char* s, ROperator op)
+int SearchOperator(char*s,ROperator op)
{
char opc;
- switch (op) {
- case ErrOp:
- case Num:
- case Var:
- return -1;
- case Juxt:
- opc = ',';
- break;
- case Add:
- opc = '+';
- break;
- case Sub:
- opc = '-';
- break;
- case Mult:
- opc = '*';
- break;
- case Div:
- opc = '/';
- break;
- case Pow:
- opc = '^';
- break;
- case NthRoot:
- opc = '#';
- break;
- case E10:
- opc = 'E';
- break;
- default:
- return -1;
+ switch(op){
+ case ErrOp:case Num:case Var:return -1;
+ case Juxt:opc=',';break;
+ case Add:opc='+';break;
+ case Sub:opc='-';break;
+ case Mult:opc='*';break;
+ case Div:opc='/';break;
+ case Pow:opc='^';break;
+ case NthRoot:opc='#';break;
+ case E10:opc='E';break;
+ default:return -1;
};
int i;
- for (i = (int)strlen(s) - 1; i >= 0; i--) {
- if (s[i] == opc && (op != Sub || (i && s[i - 1] == ')')))
- return i;
- if (s[i] == ')') {
- i = SearchCorClosebracket(s, i);
- if (i == -1)
- return -1;
- };
+ for(i=(int)strlen(s)-1;i>=0;i--){
+ if(s[i]==opc&&(op!=Sub||(i&&s[i-1]==')')))return i;
+ if(s[i]==')'){i=SearchCorClosebracket(s,i);if(i==-1)return -1;};
};
return -1;
}
-void SimplifyStr(char*& s) // Warning : deletes the old string
+void SimplifyStr(char*&s) //Warning : deletes the old string
+{if(!strlen(s))return;
+char*s1=s,*s2=s+strlen(s);signed char ind=0;
+if(s1[0]=='('&&SearchCorOpenbracket(s1,0)==s2-s1-1){
+ s1++;s2--;ind=1;}
+if(s1==s2)
{
- if (!strlen(s))
- return;
- char *s1 = s, *s2 = s + strlen(s);
- signed char ind = 0;
- if (s1[0] == '(' && SearchCorOpenbracket(s1, 0) == s2 - s1 - 1) {
- s1++;
- s2--;
- ind = 1;
- }
- if (s1 == s2) {
- delete[] s;
- s = new char[1]; // ISO C++ forbids initialization in array new
- s[0] = 0;
- return;
- }
- if (s1[0] == ' ') {
- ind = 1;
- while (s1[0] == ' ' && s1 < s2)
- s1++;
- }
- if (s1 == s2) {
- delete[] s;
- s = new char[1]; // ISO C++ forbids initialization in array new
- s[0] = 0;
- return;
- }
- if (*(s2 - 1) == ' ') {
- ind = 1;
- while (s2 > s1 && *(s2 - 1) == ' ')
- s2--;
- }
- *s2 = 0;
- s1 = CopyStr(s1);
- delete[] s;
- s = s1;
- if (ind)
- SimplifyStr(s);
+ delete[]s;
+ s=new char[1]; // ISO C++ forbids initialization in array new
+ s[0]=0;
+ return;
}
-
-int max(int a, int b)
+if(s1[0]==' '){ind=1;while(s1[0]==' '&&s1<s2)s1++;}
+if(s1==s2)
{
- return (a > b ? a : b);
+ delete[]s;
+ s=new char[1]; // ISO C++ forbids initialization in array new
+ s[0]=0;
+ return;
}
+if(*(s2-1)==' '){ind=1;while(s2>s1&&*(s2-1)==' ')s2--;}
+*s2=0;
+s1=CopyStr(s1);delete[]s;s=s1;
+if(ind)SimplifyStr(s);
+}
+
+int max(int a, int b){return (a>b?a:b);}
-int IsVar(const char* s, int n, int nvar, PRVar* ppvar)
+int IsVar(const char*s,int n,int nvar,PRVar*ppvar)
{
- if (n < 0 || n > (int)strlen(s))
- return 0;
- int i;
- int l = 0;
- for (i = 0; i < nvar; i++)
- if (CompStr(s, n, (*(ppvar + i))->name))
- l = max(l, strlen((*(ppvar + i))->name));
+ if(n<0||n>(int)strlen(s))return 0;
+ int i;int l=0;
+ for(i=0;i<nvar;i++)if(CompStr(s,n,(*(ppvar+i))->name))l=max(l,strlen((*(ppvar+i))->name));
return l;
}
-int IsFunction(const char* s, int n)
+int IsFunction(const char*s,int n)
{
- if (CompStr(s, n, "sin") || CompStr(s, n, "cos") || CompStr(s, n, "exp") || CompStr(s, n, "tan")
- || CompStr(s, n, "log")
- || CompStr(s, n, "atg")
- || CompStr(s, n, "abs"))
- return 3;
- if (CompStr(s, n, "tg") || CompStr(s, n, "ln"))
- return 2;
- if (CompStr(s, n, "sqrt") || CompStr(s, n, "asin") || CompStr(s, n, "atan") || CompStr(s, n, "acos"))
- return 4;
- if (CompStr(s, n, "arcsin") || CompStr(s, n, "arccos") || CompStr(s, n, "arctan"))
- return 6;
- if (CompStr(s, n, "arctg"))
- return 5;
+ if(CompStr(s,n,"sin")||CompStr(s,n,"cos")||CompStr(s,n,"exp")
+ ||CompStr(s,n,"tan")||CompStr(s,n,"log")||CompStr(s,n,"atg")
+ ||CompStr(s,n,"abs"))return 3;
+ if(CompStr(s,n,"tg")||CompStr(s,n,"ln"))return 2;
+ if(CompStr(s,n,"sqrt")||CompStr(s,n,"asin")||CompStr(s,n,"atan")||
+ CompStr(s,n,"acos"))return 4;
+ if(CompStr(s,n,"arcsin")||CompStr(s,n,"arccos")||CompStr(s,n,"arctan"))return 6;
+ if(CompStr(s,n,"arctg"))return 5;
return 0;
}
-int IsFunction(const char* s, int n, int nfunc, PRFunction* ppfunc)
-// Not recognized if a user-defined function is eg "sine" ie begins like
-// a standard function
-// IF PATCHED TO DO OTHERWISE, SHOULD BE PATCHED TOGETHER WITH THE
-// PARSER BELOW which treats standard functions before user-defined ones
+int IsFunction(const char*s,int n,int nfunc,PRFunction*ppfunc)
+ //Not recognized if a user-defined function is eg "sine" ie begins like
+ //a standard function
+ //IF PATCHED TO DO OTHERWISE, SHOULD BE PATCHED TOGETHER WITH THE
+ //PARSER BELOW which treats standard functions before user-defined ones
{
- int l = IsFunction(s, n);
- if (l)
- return l;
- int i;
- l = 0;
- for (i = 0; i < nfunc; i++)
- if (CompStr(s, n, ppfunc[i]->name))
- l = max(l, strlen(ppfunc[i]->name));
+ int l=IsFunction(s,n);
+ if(l)return l;
+ int i;l=0;
+ for(i=0;i<nfunc;i++)if(CompStr(s,n,ppfunc[i]->name))l=max(l,strlen(ppfunc[i]->name));
return l;
}
signed char IsFunction(ROperator op)
-{
- return (op == Exp || op == Abs || op == Sin || op == Cos || op == Tg || op == Ln || op == Atan || op == Asin
- || op == Acos
- || op == Atan
- || op == Sqrt
- || op == Opp);
+{return (op==Exp||op==Abs||op==Sin||op==Cos||op==Tg||op==Ln||
+ op==Atan||op==Asin||op==Acos||op==Atan||op==Sqrt||op==Opp);
}
-void IsolateVars(char*& s, int nvar, PRVar* ppvar, int nfunc, PRFunction* ppfunc) // Deletes the old string
+void IsolateVars(char*&s,int nvar,PRVar*ppvar,int nfunc,PRFunction*ppfunc)//Deletes the old string
{
- int i, j;
- i = 0;
- for (i = 0; s[i]; i++) {
- if (s[i] == '(') {
- i = SearchCorOpenbracket(s, i);
- if (i == -1)
- return;
- continue;
- };
- if (((j = IsVar(s, i, nvar, ppvar)) > IsFunction(s, i, nfunc, ppfunc))
- || ((CompStr(s, i, "pi") || CompStr(s, i, "PI") || CompStr(s, i, "Pi")) && (j = 2))) {
- InsStr(s, i, '(');
- InsStr(s, i + j + 1, ')');
- i += j + 1;
- continue;
- };
- if (IsFunction(s, i, nfunc, ppfunc)) {
- i += IsFunction(s, i, nfunc, ppfunc) - 1;
- if (!s[i])
- return;
- continue;
- };
+ int i,j;
+ i=0;
+ for(i=0;s[i];i++){
+ if(s[i]=='('){i=SearchCorOpenbracket(s,i);if(i==-1)return;continue;};
+ if(((j=IsVar(s,i,nvar,ppvar))>IsFunction(s,i,nfunc,ppfunc))||((CompStr(s,i,"pi")||CompStr(s,i,"PI")||CompStr(s,i,"Pi"))&&(j=2))){
+ InsStr(s,i,'(');InsStr(s,i+j+1,')');i+=j+1;continue;};
+ if(IsFunction(s,i,nfunc,ppfunc)){i+=IsFunction(s,i,nfunc,ppfunc)-1;if(!s[i])return;continue;};
};
}
-void IsolateNumbers(char*& s, int nvar, RVar** ppvar, int nfunc, RFunction** ppfunc) // Deletes the old string
+void IsolateNumbers(char*&s,int nvar,RVar**ppvar,int nfunc,RFunction**ppfunc)//Deletes the old string
{
int i, i2 = 0, ind = 0, t1, t2;
for (i = 0; s[i]; i++) {
@@ -909,10 +523,10 @@ void IsolateNumbers(char*& s, int nvar, RVar** ppvar, int nfunc, RFunction** ppf
ind = 0;
InsStr(s, i2, '(');
i++;
- InsStr(s, i, ')');
+ InsStr(s,i,')');
continue;
};
- t1 = IsVar(s, i, nvar, ppvar);
+ t1 = IsVar(s, i, nvar, ppvar );
t2 = IsFunction(s, i, nfunc, ppfunc);
if (t1 || t2) {
i += max(t1, t2) - 1;
@@ -925,387 +539,178 @@ void IsolateNumbers(char*& s, int nvar, RVar** ppvar, int nfunc, RFunction** ppf
continue;
}
if (!ind && IsNumeric(s[i])) {
- i2 = i;
+ i2 = i;
ind = 1;
};
}
- if (ind)
- InsStr(s, i2, '(');
- i++;
- InsStr(s, i, ')');
+if(ind)InsStr(s,i2,'(');i++;InsStr(s,i,')');
}
-ROperation::ROperation(const char* sp, int nvar, PRVar* ppvarp, int nfuncp, PRFunction* ppfuncp)
+ROperation::ROperation(const char*sp,int nvar,PRVar*ppvarp,int nfuncp,PRFunction*ppfuncp)
{
- ValC = ErrVal;
- mmb1 = NULL;
- mmb2 = NULL;
- pvar = NULL;
- op = ErrOp;
- pvarval = NULL;
- containfuncflag = 0;
- pfunc = NULL;
- pinstr = NULL;
- pvals = NULL;
- ppile = NULL;
- pfuncpile = NULL;
- int i, j, k, l;
- signed char flag = 1;
- char *s = CopyStr(sp), *s1 = NULL, *s2 = NULL;
- SimplifyStr(s);
- if (!s[0] || !strcmp(s, "Error")) {
- goto fin;
+ ValC=ErrVal;mmb1=NULL;mmb2=NULL;pvar=NULL;op=ErrOp;pvarval=NULL;containfuncflag=0;pfunc=NULL;pinstr=NULL;pvals=NULL;ppile=NULL;pfuncpile=NULL;
+ int i,j,k,l;signed char flag=1;
+ char*s=CopyStr(sp),*s1=NULL,*s2=NULL;
+ SimplifyStr(s);if(!s[0]||!strcmp(s,"Error")){goto fin;}
+ while(s[0]==':'||s[0]==';'){
+ s1=CopyStr(s+1);delete[]s;s=s1;s1=NULL;
+ SimplifyStr(s);if(!s[0]||!strcmp(s,"Error")){goto fin;}
}
- while (s[0] == ':' || s[0] == ';') {
- s1 = CopyStr(s + 1);
- delete[] s;
- s = s1;
- s1 = NULL;
- SimplifyStr(s);
- if (!s[0] || !strcmp(s, "Error")) {
- goto fin;
+ if(IsTNumeric(s)){op=Num;ValC=atof(s);mmb1=NULL;mmb2=NULL;goto fin;};
+ if(EqStr(s,"pi")||EqStr(s,"PI")||EqStr(s,"Pi"))
+ {op=Num;ValC=3.141592653589793238462643383279L;mmb1=NULL;mmb2=NULL;goto fin;};
+ if(IsFunction(s,0,nfuncp,ppfuncp)<IsVar(s,0,nvar,ppvarp))
+ for(i=0;i<nvar;i++)if(EqStr(s,(*(ppvarp+i))->name))
+ {pvar=ppvarp[i];pvarval=pvar->pval;op=Var;mmb1=NULL;mmb2=NULL;goto fin;};
+ for(k=0;s[k];k++){
+ if(s[k]=='('){k=SearchCorOpenbracket(s,k);if(k==-1)break;continue;};
+ if((l=IsFunction(s,k,nfuncp,ppfuncp))&&l>=IsVar(s,k,nvar,ppvarp)){
+ i=k+l;while(s[i]==' ')i++;
+ if(s[i]=='('){
+ j=SearchCorOpenbracket(s,i);
+ if(j!=-1){InsStr(s,i,';');k=j+1;}else break;
+ }else if(s[i]!=':'&&s[i]!=';'){InsStr(s,i,':');k=i;}
}
}
- if (IsTNumeric(s)) {
- op = Num;
- ValC = atof(s);
- mmb1 = NULL;
- mmb2 = NULL;
- goto fin;
- };
- if (EqStr(s, "pi") || EqStr(s, "PI") || EqStr(s, "Pi")) {
- op = Num;
- ValC = 3.141592653589793238462643383279L;
- mmb1 = NULL;
- mmb2 = NULL;
- goto fin;
- };
- if (IsFunction(s, 0, nfuncp, ppfuncp) < IsVar(s, 0, nvar, ppvarp))
- for (i = 0; i < nvar; i++)
- if (EqStr(s, (*(ppvarp + i))->name)) {
- pvar = ppvarp[i];
- pvarval = pvar->pval;
- op = Var;
- mmb1 = NULL;
- mmb2 = NULL;
- goto fin;
- };
- for (k = 0; s[k]; k++) {
- if (s[k] == '(') {
- k = SearchCorOpenbracket(s, k);
- if (k == -1)
- break;
- continue;
- };
- if ((l = IsFunction(s, k, nfuncp, ppfuncp)) && l >= IsVar(s, k, nvar, ppvarp)) {
- i = k + l;
- while (s[i] == ' ')
- i++;
- if (s[i] == '(') {
- j = SearchCorOpenbracket(s, i);
- if (j != -1) {
- InsStr(s, i, ';');
- k = j + 1;
- } else
- break;
- } else if (s[i] != ':' && s[i] != ';') {
- InsStr(s, i, ':');
- k = i;
- }
- }
- }
- IsolateNumbers(s, nvar, ppvarp, nfuncp, ppfuncp);
- if (nvar)
- IsolateVars(s, nvar, ppvarp, nfuncp, ppfuncp);
+ IsolateNumbers(s,nvar,ppvarp,nfuncp,ppfuncp);
+ if(nvar)IsolateVars(s,nvar,ppvarp,nfuncp,ppfuncp);
SupprSpaces(s);
- i = SearchOperator(s, Juxt);
- if (i != -1) {
- s1 = MidStr(s, 0, i - 1);
- s2 = MidStr(s, i + 1, strlen(s) - 1);
- op = Juxt;
- mmb1 = new ROperation(s1, nvar, ppvarp, nfuncp, ppfuncp);
- mmb2 = new ROperation(s2, nvar, ppvarp, nfuncp, ppfuncp);
- goto fin;
+ i=SearchOperator(s,Juxt);
+ if(i!=-1){s1=MidStr(s,0,i-1);s2=MidStr(s,i+1,strlen(s)-1);
+ op=Juxt;mmb1=new ROperation(s1,nvar,ppvarp,nfuncp,ppfuncp);
+ mmb2=new ROperation(s2,nvar,ppvarp,nfuncp,ppfuncp);goto fin;
};
- i = SearchOperator(s, Add);
- if (i != -1) {
- s1 = MidStr(s, 0, i - 1);
- s2 = MidStr(s, i + 1, strlen(s) - 1);
- op = Add;
- mmb1 = new ROperation(s1, nvar, ppvarp, nfuncp, ppfuncp);
- mmb2 = new ROperation(s2, nvar, ppvarp, nfuncp, ppfuncp);
- goto fin;
+ i=SearchOperator(s,Add);
+ if(i!=-1){s1=MidStr(s,0,i-1);s2=MidStr(s,i+1,strlen(s)-1);
+ op=Add;mmb1=new ROperation(s1,nvar,ppvarp,nfuncp,ppfuncp);
+ mmb2=new ROperation(s2,nvar,ppvarp,nfuncp,ppfuncp);goto fin;
};
- i = SearchOperator(s, Sub);
- if (i != -1) {
- s1 = MidStr(s, 0, i - 1);
- s2 = MidStr(s, i + 1, strlen(s) - 1);
- op = Sub;
- mmb1 = new ROperation(s1, nvar, ppvarp, nfuncp, ppfuncp);
- mmb2 = new ROperation(s2, nvar, ppvarp, nfuncp, ppfuncp);
- goto fin;
+ i=SearchOperator(s,Sub);
+ if(i!=-1){
+ s1=MidStr(s,0,i-1);s2=MidStr(s,i+1,strlen(s)-1);
+ op=Sub;mmb1=new ROperation(s1,nvar,ppvarp,nfuncp,ppfuncp);
+ mmb2=new ROperation(s2,nvar,ppvarp,nfuncp,ppfuncp);goto fin;
};
- if (s[0] == '-') {
- s2 = MidStr(s, 1, strlen(s) - 1);
- op = Opp;
- mmb1 = NULL;
- mmb2 = new ROperation(s2, nvar, ppvarp, nfuncp, ppfuncp);
- goto fin;
+ if(s[0]=='-'){s2=MidStr(s,1,strlen(s)-1);
+ op=Opp;mmb1=NULL;
+ mmb2=new ROperation(s2,nvar,ppvarp,nfuncp,ppfuncp);goto fin;
};
- for (i = 0; s[i]; i++) {
- if (s[i] == '(') {
- i = SearchCorOpenbracket(s, i);
- if (i == -1)
- break;
- continue;
- };
- if (IsFunction(s, i, nfuncp, ppfuncp)) {
- k = i + IsFunction(s, i, nfuncp, ppfuncp);
- while (s[k] == ' ')
- k++;
- if (s[k] == ';') {
+ for(i=0;s[i];i++){
+ if(s[i]=='('){i=SearchCorOpenbracket(s,i);if(i==-1)break;continue;};
+ if(IsFunction(s,i,nfuncp,ppfuncp)){
+ k=i+IsFunction(s,i,nfuncp,ppfuncp);while(s[k]==' ')k++;
+ if(s[k]==';'){
// s=DelStr(s,k);
- j = k;
- while (s[j] != '(')
- j++;
- j = SearchCorOpenbracket(s, j);
- if (j != -1) {
- InsStr(s, j, ')');
- InsStr(s, i, '(');
- i = j + 2;
- }
- } else if (s[k] == ':') {
+ j=k;while(s[j]!='(')j++;
+ j=SearchCorOpenbracket(s,j);
+ if(j!=-1){InsStr(s,j,')');InsStr(s,i,'(');i=j+2;}
+ }else if(s[k]==':'){
// s=DelStr(s,k);
- for (j = k; s[j]; j++)
- if (s[j] == '(') {
- j = SearchCorOpenbracket(s, j);
- break;
- }
- if (j == -1)
- break;
- for (j++; s[j]; j++) {
- if (s[j] == '(') {
- j = SearchCorOpenbracket(s, j);
- if (j == -1) {
- flag = 0;
- break;
- };
- continue;
- };
- if (IsFunction(s, j, nfuncp, ppfuncp))
- break;
- }
- if (flag == 0) {
- flag = 1;
- break;
+ for(j=k;s[j];j++)
+ if(s[j]=='('){j=SearchCorOpenbracket(s,j);break;}
+ if(j==-1)break;
+ for(j++;s[j];j++){
+ if(s[j]=='('){j=SearchCorOpenbracket(s,j);if(j==-1){flag=0;break;};continue;};
+ if(IsFunction(s,j,nfuncp,ppfuncp))break;
}
- while (j > i && s[j - 1] != ')')
- j--;
- if (j <= i + 1)
- break;
- InsStr(s, i, '(');
- InsStr(s, j + 1, ')');
- i = j + 1;
+ if(flag==0){flag=1;break;}
+ while(j>i&&s[j-1]!=')')j--;if(j<=i+1)break;
+ InsStr(s,i,'(');InsStr(s,j+1,')');
+ i=j+1;
}
}
}
- for (i = 0; s[i] && s[i + 1]; i++)
- if (s[i] == ')' && s[i + 1] == '(')
- InsStr(s, ++i, '*');
- if (s[0] == '(' && SearchCorOpenbracket(s, 0) == (int)strlen(s) - 1) {
- if (CompStr(s, 1, "exp")) {
- op = Exp;
- s2 = MidStr(s, 4, strlen(s) - 2);
- } else if (CompStr(s, 1, "abs")) {
- op = Abs;
- s2 = MidStr(s, 4, strlen(s) - 2);
- } else if (CompStr(s, 1, "sin")) {
- op = Sin;
- s2 = MidStr(s, 4, strlen(s) - 2);
- } else if (CompStr(s, 1, "cos")) {
- op = Cos;
- s2 = MidStr(s, 4, strlen(s) - 2);
- } else if (CompStr(s, 1, "tan")) {
- op = Tg;
- s2 = MidStr(s, 4, strlen(s) - 2);
- } else if (CompStr(s, 1, "log")) {
- op = Ln;
- s2 = MidStr(s, 4, strlen(s) - 2);
- } else if (CompStr(s, 1, "atg")) {
- op = Atan;
- s2 = MidStr(s, 4, strlen(s) - 2);
- } else if (CompStr(s, 1, "tg")) {
- op = Tg;
- s2 = MidStr(s, 3, strlen(s) - 2);
- } else if (CompStr(s, 1, "ln")) {
- op = Ln;
- s2 = MidStr(s, 3, strlen(s) - 2);
- } else if (CompStr(s, 1, "asin")) {
- op = Asin;
- s2 = MidStr(s, 5, strlen(s) - 2);
- } else if (CompStr(s, 1, "acos")) {
- op = Acos;
- s2 = MidStr(s, 5, strlen(s) - 2);
- } else if (CompStr(s, 1, "atan")) {
- op = Atan;
- s2 = MidStr(s, 5, strlen(s) - 2);
- } else if (CompStr(s, 1, "sqrt")) {
- op = Sqrt;
- s2 = MidStr(s, 5, strlen(s) - 2);
- } else if (CompStr(s, 1, "arcsin")) {
- op = Asin;
- s2 = MidStr(s, 7, strlen(s) - 2);
- } else if (CompStr(s, 1, "arccos")) {
- op = Acos;
- s2 = MidStr(s, 7, strlen(s) - 2);
- } else if (CompStr(s, 1, "arctan")) {
- op = Atan;
- s2 = MidStr(s, 7, strlen(s) - 2);
- } else if (CompStr(s, 1, "arctg")) {
- op = Atan;
- s2 = MidStr(s, 6, strlen(s) - 2);
- } else {
- for (i = -1, k = 0, j = 0; j < nfuncp; j++)
- if (CompStr(s, 1, ppfuncp[j]->name) && k < (int)strlen(ppfuncp[j]->name)) {
- k = strlen(ppfuncp[j]->name);
- i = j;
- }
- if (i > -1) {
- op = Fun;
- s2 = MidStr(s, strlen(ppfuncp[i]->name) + 1, strlen(s) - 2);
- pfunc = ppfuncp[i];
+ for(i=0;s[i]&&s[i+1];i++)if(s[i]==')'&&s[i+1]=='(')
+ InsStr(s,++i,'*');
+ if(s[0]=='('&&SearchCorOpenbracket(s,0)==(int)strlen(s)-1){
+ if(CompStr(s,1,"exp")){op=Exp;s2=MidStr(s,4,strlen(s)-2);}
+ else if(CompStr(s,1,"abs")){op=Abs;s2=MidStr(s,4,strlen(s)-2);}
+ else if(CompStr(s,1,"sin")){op=Sin;s2=MidStr(s,4,strlen(s)-2);}
+ else if(CompStr(s,1,"cos")){op=Cos;s2=MidStr(s,4,strlen(s)-2);}
+ else if(CompStr(s,1,"tan")){op=Tg;s2=MidStr(s,4,strlen(s)-2);}
+ else if(CompStr(s,1,"log")){op=Ln;s2=MidStr(s,4,strlen(s)-2);}
+ else if(CompStr(s,1,"atg")){op=Atan;s2=MidStr(s,4,strlen(s)-2);}
+ else if(CompStr(s,1,"tg")){op=Tg;s2=MidStr(s,3,strlen(s)-2);}
+ else if(CompStr(s,1,"ln")){op=Ln;s2=MidStr(s,3,strlen(s)-2);}
+ else if(CompStr(s,1,"asin")){op=Asin;s2=MidStr(s,5,strlen(s)-2);}
+ else if(CompStr(s,1,"acos")){op=Acos;s2=MidStr(s,5,strlen(s)-2);}
+ else if(CompStr(s,1,"atan")){op=Atan;s2=MidStr(s,5,strlen(s)-2);}
+ else if(CompStr(s,1,"sqrt")){op=Sqrt;s2=MidStr(s,5,strlen(s)-2);}
+ else if(CompStr(s,1,"arcsin")){op=Asin;s2=MidStr(s,7,strlen(s)-2);}
+ else if(CompStr(s,1,"arccos")){op=Acos;s2=MidStr(s,7,strlen(s)-2);}
+ else if(CompStr(s,1,"arctan")){op=Atan;s2=MidStr(s,7,strlen(s)-2);}
+ else if(CompStr(s,1,"arctg")){op=Atan;s2=MidStr(s,6,strlen(s)-2);}
+ else {
+ for(i=-1,k=0,j=0;j<nfuncp;j++)if(CompStr(s,1,ppfuncp[j]->name)&&k<(int)strlen(ppfuncp[j]->name)){k=strlen(ppfuncp[j]->name);i=j;}
+ if(i>-1){
+ op=Fun;s2=MidStr(s,strlen(ppfuncp[i]->name)+1,strlen(s)-2);
+ pfunc=ppfuncp[i];
}
}
- mmb1 = NULL;
- mmb2 = new ROperation(s2, nvar, ppvarp, nfuncp, ppfuncp);
- if (op == Fun)
- if (mmb2->NMembers() != pfunc->nvars) {
- op = ErrOp;
- mmb1 = NULL;
- mmb2 = NULL;
- goto fin;
- }
+ mmb1=NULL;mmb2=new ROperation(s2,nvar,ppvarp,nfuncp,ppfuncp);
+ if(op==Fun)if(mmb2->NMembers()!=pfunc->nvars){op=ErrOp;mmb1=NULL;mmb2=NULL;goto fin;}
goto fin;
};
- i = SearchOperator(s, Mult);
- if (i != -1) {
- s1 = MidStr(s, 0, i - 1);
- s2 = MidStr(s, i + 1, strlen(s) - 1);
- op = Mult;
- mmb1 = new ROperation(s1, nvar, ppvarp, nfuncp, ppfuncp);
- mmb2 = new ROperation(s2, nvar, ppvarp, nfuncp, ppfuncp);
- goto fin;
+ i=SearchOperator(s,Mult);
+ if(i!=-1){s1=MidStr(s,0,i-1);s2=MidStr(s,i+1,strlen(s)-1);
+ op=Mult;mmb1=new ROperation(s1,nvar,ppvarp,nfuncp,ppfuncp);
+ mmb2=new ROperation(s2,nvar,ppvarp,nfuncp,ppfuncp);goto fin;
};
- i = SearchOperator(s, Div);
- if (i != -1) {
- s1 = MidStr(s, 0, i - 1);
- s2 = MidStr(s, i + 1, strlen(s) - 1);
- op = Div;
- mmb1 = new ROperation(s1, nvar, ppvarp, nfuncp, ppfuncp);
- mmb2 = new ROperation(s2, nvar, ppvarp, nfuncp, ppfuncp);
- goto fin;
+ i=SearchOperator(s,Div);
+ if(i!=-1){s1=MidStr(s,0,i-1);s2=MidStr(s,i+1,strlen(s)-1);
+ op=Div;mmb1=new ROperation(s1,nvar,ppvarp,nfuncp,ppfuncp);
+ mmb2=new ROperation(s2,nvar,ppvarp,nfuncp,ppfuncp);goto fin;
};
- i = SearchOperator(s, Pow);
- if (i != -1) {
- s1 = MidStr(s, 0, i - 1);
- s2 = MidStr(s, i + 1, strlen(s) - 1);
- op = Pow;
- mmb1 = new ROperation(s1, nvar, ppvarp, nfuncp, ppfuncp);
- mmb2 = new ROperation(s2, nvar, ppvarp, nfuncp, ppfuncp);
- goto fin;
+ i=SearchOperator(s,Pow);
+ if(i!=-1){s1=MidStr(s,0,i-1);s2=MidStr(s,i+1,strlen(s)-1);
+ op=Pow;mmb1=new ROperation(s1,nvar,ppvarp,nfuncp,ppfuncp);
+ mmb2=new ROperation(s2,nvar,ppvarp,nfuncp,ppfuncp);goto fin;
};
- i = SearchOperator(s, NthRoot);
- if (i != -1) {
- s1 = MidStr(s, 0, i - 1);
- s2 = MidStr(s, i + 1, strlen(s) - 1);
- if (i == 0 || s[i - 1] != ')') {
- op = Sqrt;
- mmb1 = NULL;
- } else {
- op = NthRoot;
- mmb1 = new ROperation(s1, nvar, ppvarp, nfuncp, ppfuncp);
- };
- mmb2 = new ROperation(s2, nvar, ppvarp, nfuncp, ppfuncp);
- goto fin;
+ i=SearchOperator(s,NthRoot);
+ if(i!=-1){s1=MidStr(s,0,i-1);s2=MidStr(s,i+1,strlen(s)-1);
+ if(i==0||s[i-1]!=')')
+ {op=Sqrt;mmb1=NULL;}else
+ {op=NthRoot;mmb1=new ROperation(s1,nvar,ppvarp,nfuncp,ppfuncp);};
+ mmb2=new ROperation(s2,nvar,ppvarp,nfuncp,ppfuncp);goto fin;
};
- i = SearchOperator(s, E10);
- if (i != -1) {
- s1 = MidStr(s, 0, i - 1);
- s2 = MidStr(s, i + 1, strlen(s) - 1);
- op = E10;
- mmb1 = new ROperation(s1, nvar, ppvarp, nfuncp, ppfuncp);
- mmb2 = new ROperation(s2, nvar, ppvarp, nfuncp, ppfuncp);
- goto fin;
+ i=SearchOperator(s,E10);
+ if(i!=-1){s1=MidStr(s,0,i-1);s2=MidStr(s,i+1,strlen(s)-1);
+ op=E10;mmb1=new ROperation(s1,nvar,ppvarp,nfuncp,ppfuncp);
+ mmb2=new ROperation(s2,nvar,ppvarp,nfuncp,ppfuncp);goto fin;
};
- op = ErrOp;
- mmb1 = NULL;
- mmb2 = NULL;
+ op=ErrOp;mmb1=NULL;mmb2=NULL;
fin:
BuildCode();
- delete[] s;
- if (s1 != NULL)
- delete[] s1;
- if (s2 != NULL)
- delete[] s2;
+ delete[]s;if(s1!=NULL)delete[] s1;if(s2!=NULL)delete[]s2;
}
void ROperation::Destroy()
{
- if (mmb1 != NULL && mmb2 != NULL && mmb1 != mmb2) {
- delete mmb1;
- delete mmb2;
- mmb1 = NULL;
- mmb2 = NULL;
- } else if (mmb1 != NULL) {
- delete mmb1;
- mmb1 = NULL;
- } else if (mmb2 != NULL) {
- delete mmb2;
- mmb2 = NULL;
- }
- if (pinstr != NULL) {
- delete[] pinstr;
- pinstr = NULL;
- }
- if (pvals != NULL) {
- if (op == ErrOp || op == Num)
- delete pvals[0];
- delete[] pvals;
- pvals = NULL;
- }
- if (ppile != NULL) {
- delete[] ppile;
- ppile = NULL;
- }
- if (pfuncpile != NULL) {
- delete[] pfuncpile;
- pfuncpile = NULL;
+ if(mmb1!=NULL&&mmb2!=NULL&&mmb1!=mmb2){delete mmb1;delete mmb2;mmb1=NULL;mmb2=NULL;}
+ else if(mmb1!=NULL){delete mmb1;mmb1=NULL;}else if(mmb2!=NULL){delete mmb2;mmb2=NULL;}
+ if(pinstr!=NULL){delete[]pinstr;pinstr=NULL;}
+ if(pvals!=NULL){
+ if(op==ErrOp||op==Num)delete pvals[0];
+ delete[]pvals;pvals=NULL;
}
+ if(ppile!=NULL){delete[]ppile;ppile=NULL;}
+ if(pfuncpile!=NULL){delete[]pfuncpile;pfuncpile=NULL;}
}
-int operator==(const RVar& var1, const RVar& var2)
-{
- return (var1.pval == var2.pval && EqStr(var1.name, var2.name));
+int operator==(const RVar& var1,const RVar& var2)
+{return(var1.pval==var2.pval&&EqStr(var1.name,var2.name));
}
-int operator==(const RFunction& f1, const RFunction& f2)
+int operator==(const RFunction& f1,const RFunction& f2)
{
- if (f1.type != f2.type)
- return 0;
- if (f1.type == -1)
- return 1; // Nonfunction==nonfunction
- if (f1.type == 0)
- return (f1.pfuncval == f2.pfuncval && EqStr(f1.name, f2.name));
- if (f1.op != f2.op)
- return 0;
- if (!EqStr(f1.name, f2.name))
- return 0;
- if (f1.nvars != f2.nvars)
- return 0;
+ if(f1.type!=f2.type)return 0;
+ if(f1.type==-1)return 1; // Nonfunction==nonfunction
+ if(f1.type==0)return (f1.pfuncval==f2.pfuncval&&EqStr(f1.name,f2.name));
+ if(f1.op!=f2.op)return 0;
+ if(!EqStr(f1.name,f2.name))return 0;
+ if(f1.nvars!=f2.nvars)return 0;
int i;
- for (i = 0; i < f1.nvars; i++)
- if (!(*f1.ppvar[i] == *f2.ppvar[i]))
- return 0;
+ for(i=0;i<f1.nvars;i++)if(!(*f1.ppvar[i]==*f2.ppvar[i]))return 0;
return 1;
}
@@ -1313,10 +718,8 @@ int operator==(const RFunction& f1, const RFunction& f2)
double ROperation::Val() const // Won't work if multi-variable functions are included
{
double v1=ErrVal,v2=ErrVal;
- if(mmb1!=NULL){v1=mmb1->Val();if(fabsl(v1)<sqrtminfloat)v1=0;else if(v1==ErrVal||fabsl(v1)>sqrtmaxfloat)return
-ErrVal;};
- if(mmb2!=NULL){v2=mmb2->Val();if(fabsl(v2)<sqrtminfloat)v2=0;else if(v2==ErrVal||fabsl(v2)>sqrtmaxfloat)return
-ErrVal;};
+ if(mmb1!=NULL){v1=mmb1->Val();if(fabsl(v1)<sqrtminfloat)v1=0;else if(v1==ErrVal||fabsl(v1)>sqrtmaxfloat)return ErrVal;};
+ if(mmb2!=NULL){v2=mmb2->Val();if(fabsl(v2)<sqrtminfloat)v2=0;else if(v2==ErrVal||fabsl(v2)>sqrtmaxfloat)return ErrVal;};
switch(op){
case Num:return ValC;
case Var:return *pvarval;
@@ -1325,8 +728,7 @@ ErrVal;};
case Opp:return -v2;
case Mult:return v1*v2;
case Div:if(v2)return v1/v2;else return ErrVal;
- case Pow:if(v1==0)return 0;else if((v1>0||!fmodl(v2,1))&&v2*logl(fabsl(v1))<DBL_MAX_EXP)return powl(v1,v2);else return
-ErrVal;
+ case Pow:if(v1==0)return 0;else if((v1>0||!fmodl(v2,1))&&v2*logl(fabsl(v1))<DBL_MAX_EXP)return powl(v1,v2);else return ErrVal;
case Sqrt:if(v2>=0)return sqrtl(v2);else return ErrVal;
case NthRoot:if(!v1||v2*logl(fabsl(v1))<DBL_MIN_EXP)return ErrVal;
else if(v2>=0)return powl(v2,1/v1);else
@@ -1338,8 +740,7 @@ ErrVal;
case Cos:if(fabsl(v2)<inveps)return cosl(v2);else return ErrVal;
case Tg:if(fabsl(v2)<inveps)return tanl(v2);else return ErrVal;
case Atan:
- if(mmb2->op==Juxt){v1=mmb2->NthMember(1).Val();v2=mmb2->NthMember(2).Val();return (v1||v2?atan2(v1,v2):ErrVal);}else
-return atanl(v2);
+ if(mmb2->op==Juxt){v1=mmb2->NthMember(1).Val();v2=mmb2->NthMember(2).Val();return (v1||v2?atan2(v1,v2):ErrVal);}else return atanl(v2);
case Asin:if(v2<-1||v2>1)return ErrVal;else return asinl(v2);
case Acos:if(v2<-1||v2>1)return ErrVal;else return acosl(v2);
case Abs:return fabsl(v2);
@@ -1350,992 +751,511 @@ return atanl(v2);
*/
signed char ROperation::ContainVar(const RVar& varp) const
-{
- if (op == Var) {
- if (EqStr(pvar->name, varp.name) && pvar->pval == varp.pval)
- return 1;
- else
- return 0;
- };
- if (mmb1 != NULL && mmb1->ContainVar(varp))
- return 1;
- if (mmb2 != NULL && mmb2->ContainVar(varp))
- return 1;
- return 0;
+{if(op==Var){if(EqStr(pvar->name,varp.name)&&pvar->pval==varp.pval)
+ return 1;else return 0;};
+if(mmb1!=NULL&&mmb1->ContainVar(varp))return 1;
+if(mmb2!=NULL&&mmb2->ContainVar(varp))return 1;
+return 0;
}
signed char ROperation::ContainFuncNoRec(const RFunction& func) const // No recursive test on subfunctions
-{
- if (op == Fun) {
- if (*pfunc == func)
- return 1;
- else
- return 0;
- }
- if (mmb1 != NULL && mmb1->ContainFuncNoRec(func))
- return 1;
- if (mmb2 != NULL && mmb2->ContainFuncNoRec(func))
- return 1;
- return 0;
+{if(op==Fun){if(*pfunc==func)
+ return 1;else return 0;}
+ if(mmb1!=NULL&&mmb1->ContainFuncNoRec(func))return 1;
+ if(mmb2!=NULL&&mmb2->ContainFuncNoRec(func))return 1;
+ return 0;
}
signed char ROperation::ContainFunc(const RFunction& func) const // Recursive test on subfunctions
{
- if (containfuncflag)
- return 0;
- if (op == Fun && *pfunc == func)
- return 1;
- containfuncflag = 1;
- if (op == Fun)
- if (pfunc->op.ContainFunc(func)) {
- containfuncflag = 0;
- return 1;
- }
- if (mmb1 != NULL && mmb1->ContainFunc(func)) {
- containfuncflag = 0;
- return 1;
- }
- if (mmb2 != NULL && mmb2->ContainFunc(func)) {
- containfuncflag = 0;
- return 1;
- }
- containfuncflag = 0;
- return 0;
+ if(containfuncflag)return 0;
+ if(op==Fun&&*pfunc==func)return 1;
+ containfuncflag=1;
+ if(op==Fun)if(pfunc->op.ContainFunc(func)){containfuncflag=0;return 1;}
+ if(mmb1!=NULL&&mmb1->ContainFunc(func)){containfuncflag=0;return 1;}
+ if(mmb2!=NULL&&mmb2->ContainFunc(func)){containfuncflag=0;return 1;}
+ containfuncflag=0;return 0;
}
-signed char ROperation::HasError(const ROperation* pop) const
+signed char ROperation::HasError(const ROperation*pop) const
{
- if (op == ErrOp)
- return 1;
- if (op == Fun && pfunc->type == 1 && pfunc->op == *(pop == NULL ? this : pop))
- return 1;
- if (op == Fun && pfunc->type == 1 && pfunc->op.HasError((pop == NULL ? this : pop)))
- return 1;
- if (mmb1 != NULL && mmb1->HasError((pop == NULL ? this : pop)))
- return 1;
- if (mmb2 != NULL && mmb2->HasError((pop == NULL ? this : pop)))
- return 1;
- if (op == Fun && pfunc->type == -1)
- return 1;
+ if(op==ErrOp)return 1;
+ if(op==Fun&&pfunc->type==1&&pfunc->op==*(pop==NULL?this:pop))return 1;
+ if(op==Fun&&pfunc->type==1&&pfunc->op.HasError((pop==NULL?this:pop)))return 1;
+ if(mmb1!=NULL&&mmb1->HasError((pop==NULL?this:pop)))return 1;
+ if(mmb2!=NULL&&mmb2->HasError((pop==NULL?this:pop)))return 1;
+ if(op==Fun&&pfunc->type==-1)return 1;
return 0;
}
-int ROperation::NMembers() const // Number of members for an operation like a,b,c...
+int ROperation::NMembers() const //Number of members for an operation like a,b,c...
{
- if (op == Fun)
- return (pfunc->type == 1 ? pfunc->op.NMembers() : pfunc->type == 0 ? 1 : 0);
- if (op != Juxt)
- return 1;
- else if (mmb2 == NULL)
- return 0;
- else
- return 1 + mmb2->NMembers();
+ if(op==Fun)return(pfunc->type==1?pfunc->op.NMembers():pfunc->type==0?1:0);
+ if(op!=Juxt)return 1;else if(mmb2==NULL)return 0;else return 1+mmb2->NMembers();
}
ROperation ROperation::NthMember(int n) const
{
PRFunction prf;
- if (op == Fun && pfunc->type == 1 && pfunc->op.NMembers() > 1) {
- prf = new RFunction(pfunc->op.NthMember(n), pfunc->nvars, pfunc->ppvar);
- char* s = new char[strlen(pfunc->name) + 10];
- sprintf(s, "(%s_%i)", pfunc->name, n);
- prf->SetName(s);
- delete[] s;
- return (*prf)(*mmb2);
+ if(op==Fun&&pfunc->type==1&&pfunc->op.NMembers()>1){
+ prf=new RFunction(pfunc->op.NthMember(n),pfunc->nvars,pfunc->ppvar);
+ char*s=new char[strlen(pfunc->name)+10];
+ sprintf(s,"(%s_%i)",pfunc->name,n);prf->SetName(s);delete[]s;
+ return(*prf)(*mmb2);
}
- if (n == 1) {
- if (op != Juxt)
- return *this;
- else if (mmb1 != NULL)
- return *mmb1;
- else
- return ErrVal;
+ if(n==1){
+ if(op!=Juxt)return *this; else if(mmb1!=NULL)return *mmb1;else return ErrVal;
};
- if (op != Juxt)
- return ErrVal;
- if (n > 1 && mmb2 != NULL)
- return mmb2->NthMember(n - 1);
+ if(op!=Juxt)return ErrVal;
+ if(n>1&&mmb2!=NULL)return mmb2->NthMember(n-1);
return ErrVal;
}
-ROperation ROperation::Substitute(const RVar& var,
- const ROperation& rop) const // Replaces variable var with expression rop
+ROperation ROperation::Substitute(const RVar& var,const ROperation& rop) const // Replaces variable var with expression rop
{
- if (!ContainVar(var))
- return *this;
- if (op == Var)
- return rop;
+ if(!ContainVar(var))return *this;
+ if(op==Var)return rop;
ROperation r;
- r.op = op;
- r.pvar = pvar;
- r.pvarval = pvarval;
- r.ValC = ValC;
- r.pfunc = pfunc;
- if (mmb1 != NULL)
- r.mmb1 = new ROperation(mmb1->Substitute(var, rop));
- else
- r.mmb1 = NULL;
- if (mmb2 != NULL)
- r.mmb2 = new ROperation(mmb2->Substitute(var, rop));
- else
- r.mmb2 = NULL;
+ r.op=op;r.pvar=pvar;r.pvarval=pvarval;r.ValC=ValC;r.pfunc=pfunc;
+ if(mmb1!=NULL)r.mmb1=new ROperation(mmb1->Substitute(var,rop));else r.mmb1=NULL;
+ if(mmb2!=NULL)r.mmb2=new ROperation(mmb2->Substitute(var,rop));else r.mmb2=NULL;
return r;
}
ROperation ROperation::Diff(const RVar& var) const
{
- if (!ContainVar(var))
- return 0.0;
- if (op == Var)
- return 1.0;
- ROperation **ppop1, op2;
- int i, j;
- switch (op) {
- case Juxt:
- return (mmb1->Diff(var), mmb2->Diff(var));
- case Add:
- return (mmb1->Diff(var) + mmb2->Diff(var));
- case Sub:
- return (mmb1->Diff(var) - mmb2->Diff(var));
- case Opp:
- return (-mmb2->Diff(var));
- case Mult:
- return ((*mmb1) * (mmb2->Diff(var)) + (*mmb2) * (mmb1->Diff(var)));
- case Div:
- if (mmb2->ContainVar(var))
- return (((*mmb2) * (mmb1->Diff(var)) - (*mmb1) * (mmb2->Diff(var))) / ((*mmb2) ^ 2));
- else
- return (mmb1->Diff(var) / (*mmb2));
- case Pow:
- if (mmb2->ContainVar(var))
- return ((*this) * (log(*mmb1) * mmb2->Diff(var) + (*mmb2) * mmb1->Diff(var) / (*mmb1)));
- else
- return (*mmb2) * mmb1->Diff(var) * ((*mmb1) ^ (*mmb2 - 1));
- case Sqrt:
- return (mmb2->Diff(var) / (2 * sqrt(*mmb2)));
- case NthRoot: {
- ROperation interm = (*mmb2) ^ (1 / (*mmb1));
- return interm.Diff(var);
- };
- case E10: {
- ROperation interm = (*mmb1) * (10 ^ (*mmb2));
- return interm.Diff(var);
- };
- ;
- case Ln:
- return (mmb2->Diff(var) / (*mmb2));
- case Exp:
- return (mmb2->Diff(var) * (*this));
- case Sin:
- return (mmb2->Diff(var) * cos(*mmb2));
- case Cos:
- return (-mmb2->Diff(var) * sin(*mmb2));
- case Tg:
- return (mmb2->Diff(var) * (1 + ((*this) ^ 2)));
- case Atan:
- if (mmb2->op != Juxt)
- return (mmb2->Diff(var) / (1 + ((*mmb2) ^ 2)));
- else
- return ((mmb2->NthMember(1).Diff(var)) * (mmb2->NthMember(2))
- - (mmb2->NthMember(2).Diff(var)) * (mmb2->NthMember(1)))
- / (((mmb2->NthMember(1)) ^ 2) + ((mmb2->NthMember(2)) ^ 2));
- case Asin:
- return (mmb2->Diff(var) / sqrt(1 - ((*mmb2) ^ 2)));
- case Acos:
- return (-mmb2->Diff(var) / sqrt(1 - ((*mmb2) ^ 2)));
- case Abs:
- return (mmb2->Diff(var) * (*mmb2) / (*this));
- case Fun:
- if (pfunc->type == -1 || pfunc->type == 0)
- return ErrVal;
- if (pfunc->nvars == 0)
- return 0.;
- else if (pfunc->op.NMembers() > 1) {
- j = pfunc->op.NMembers();
- ppop1 = new ROperation*[j];
- for (i = 0; i < j; i++)
- ppop1[i] = new ROperation(NthMember(i + 1).Diff(var));
- op2 = ApplyOperator(pfunc->nvars, ppop1, &operator, );
- for (i = 0; i < pfunc->nvars; i++)
- delete ppop1[i];
- delete[] ppop1;
- return op2;
- } else {
- ppop1 = new ROperation*[pfunc->nvars];
- for (i = 0; i < pfunc->nvars; i++) {
- ppop1[i] = new ROperation(pfunc->op.Diff(*pfunc->ppvar[i]));
- for (j = 0; j < pfunc->nvars; j++)
- *ppop1[i] = ppop1[i]->Substitute(*pfunc->ppvar[j], mmb2->NthMember(j + 1));
- *ppop1[i] = (mmb2->NthMember(i + 1).Diff(var)) * (*ppop1[i]);
- }
- op2 = ApplyOperator(pfunc->nvars, ppop1, &::operator+);
- for (i = 0; i < pfunc->nvars; i++)
- delete ppop1[i];
- delete[] ppop1;
- return op2;
- // In the obtained expression, f' will have been replaced with its expression but f will remain pointing to
- // itself ; this could cause some trouble if changing f afterwards
+ if(!ContainVar(var))return 0.0;
+ if(op==Var)return 1.0;
+ ROperation **ppop1,op2;int i,j;
+ switch(op){
+ case Juxt:return(mmb1->Diff(var),mmb2->Diff(var));
+ case Add:return(mmb1->Diff(var)+mmb2->Diff(var));
+ case Sub:return(mmb1->Diff(var)-mmb2->Diff(var));
+ case Opp:return(-mmb2->Diff(var));
+ case Mult:return((*mmb1)*(mmb2->Diff(var))+(*mmb2)*(mmb1->Diff(var)));
+ case Div:if(mmb2->ContainVar(var))return(((*mmb2)*(mmb1->Diff(var))-(*mmb1)*(mmb2->Diff(var)))/((*mmb2)^2));
+ else return(mmb1->Diff(var)/(*mmb2));
+ case Pow:if(mmb2->ContainVar(var))return((*this)*(log(*mmb1)*mmb2->Diff(var)+
+ (*mmb2)*mmb1->Diff(var)/(*mmb1)));else
+ return (*mmb2)*mmb1->Diff(var)*((*mmb1)^(*mmb2-1));
+ case Sqrt:return(mmb2->Diff(var)/(2*sqrt(*mmb2)));
+ case NthRoot:{ROperation interm=(*mmb2)^(1/(*mmb1));return interm.Diff(var);};
+ case E10:{ROperation interm=(*mmb1)*(10^(*mmb2));return interm.Diff(var);};;
+ case Ln:return (mmb2->Diff(var)/(*mmb2));
+ case Exp:return (mmb2->Diff(var)*(*this));
+ case Sin:return (mmb2->Diff(var)*cos(*mmb2));
+ case Cos:return (-mmb2->Diff(var)*sin(*mmb2));
+ case Tg:return (mmb2->Diff(var)*(1+((*this)^2)));
+ case Atan:
+ if(mmb2->op!=Juxt)return(mmb2->Diff(var)/(1+((*mmb2)^2)));
+ else return ((mmb2->NthMember(1).Diff(var))*(mmb2->NthMember(2))-(mmb2->NthMember(2).Diff(var))*(mmb2->NthMember(1)))/(((mmb2->NthMember(1))^2)+((mmb2->NthMember(2))^2));
+ case Asin:return(mmb2->Diff(var)/sqrt(1-((*mmb2)^2)));
+ case Acos:return(-mmb2->Diff(var)/sqrt(1-((*mmb2)^2)));
+ case Abs:return(mmb2->Diff(var)*(*mmb2)/(*this));
+ case Fun:if(pfunc->type==-1||pfunc->type==0)return ErrVal;
+ if(pfunc->nvars==0)return 0.;
+ else if(pfunc->op.NMembers()>1){
+ j=pfunc->op.NMembers();
+ ppop1=new ROperation*[j];
+ for(i=0;i<j;i++)ppop1[i]=new ROperation(NthMember(i+1).Diff(var));
+ op2=ApplyOperator(pfunc->nvars,ppop1,&operator,);
+ for(i=0;i<pfunc->nvars;i++)delete ppop1[i];
+ delete[]ppop1;
+ return op2;
+ }else{
+ ppop1=new ROperation*[pfunc->nvars];
+ for(i=0;i<pfunc->nvars;i++){
+ ppop1[i]=new ROperation(pfunc->op.Diff(*pfunc->ppvar[i]));
+ for(j=0;j<pfunc->nvars;j++)
+ *ppop1[i]=ppop1[i]->Substitute(*pfunc->ppvar[j],mmb2->NthMember(j+1));
+ *ppop1[i]=(mmb2->NthMember(i+1).Diff(var))*(*ppop1[i]);
}
- default:
- return ErrVal;
+ op2=ApplyOperator(pfunc->nvars,ppop1,&::operator+);
+ for(i=0;i<pfunc->nvars;i++)delete ppop1[i];
+ delete[]ppop1;
+ return op2;
+ // In the obtained expression, f' will have been replaced with its expression but f will remain pointing to itself ; this could cause some trouble if changing f afterwards
+ }
+ default:return ErrVal;
};
}
char* ValToStr(double x)
{
- char* s = new char[30];
- if (x == (double)3.141592653589793238462643383279L)
- sprintf(s, "pi");
- else
- sprintf(s, "%.16G", x);
+ char*s=new char[30];
+ if(x==(double)3.141592653589793238462643383279L)sprintf(s,"pi");else sprintf(s,"%.16G",x);
return s;
}
char* ROperation::Expr() const
{
- char *s = NULL, *s1 = NULL, *s2 = NULL;
- int n = 10;
- signed char f = 0, g = 0;
- if (op == Fun)
- if (strlen(pfunc->name) > 4)
- n += strlen(pfunc->name) - 4;
- if (mmb1 != NULL) {
- s1 = mmb1->Expr();
- n += strlen(s1);
- f = IsFunction(mmb1->op);
- }
- if (mmb2 != NULL) {
- s2 = mmb2->Expr();
- n += strlen(s2);
- g = IsFunction(mmb2->op);
- }
- s = new char[n];
- switch (op) {
- case Num:
- return ValToStr(ValC);
- case Var:
- return CopyStr(pvar->name);
- case Juxt:
- sprintf(s, "%s , %s", s1, s2);
- break;
- case Add:
- f = f || (mmb1->op == Juxt);
- g = g || (mmb2->op == Juxt);
- if (f && g)
- sprintf(s, "(%s)+(%s)", s1, s2);
- else if (f)
- sprintf(s, "(%s)+%s", s1, s2);
- else if (g)
- sprintf(s, "%s+(%s)", s1, s2);
- else
- sprintf(s, "%s+%s", s1, s2);
- break;
- case Sub:
- f = f || (mmb1->op == Juxt);
- g = g || (mmb2->op == Juxt || mmb2->op == Add || mmb2->op == Sub);
- if (f && g)
- sprintf(s, "(%s)-(%s)", s1, s2);
- else if (f)
- sprintf(s, "(%s)-%s", s1, s2);
- else if (g)
- sprintf(s, "%s-(%s)", s1, s2);
- else
- sprintf(s, "%s-%s", s1, s2);
- break;
- case Opp:
- if (mmb2->op == Add || mmb2->op == Sub || mmb2->op == Juxt)
- sprintf(s, "-(%s)", s2);
- else
- sprintf(s, "-%s", s2);
- break;
- case Mult:
- f = f || (mmb1->op == Juxt || mmb1->op == Add || mmb1->op == Sub || mmb1->op == Opp || mmb1->op == Div);
- g = g || (mmb2->op == Juxt || mmb2->op == Add || mmb2->op == Sub || mmb2->op == Opp);
- if (f && g)
- sprintf(s, "(%s)*(%s)", s1, s2);
- else if (f)
- sprintf(s, "(%s)*%s", s1, s2);
- else if (g)
- sprintf(s, "%s*(%s)", s1, s2);
- else
- sprintf(s, "%s*%s", s1, s2);
- break;
- case Div:
- f = f || (mmb1->op == Juxt || mmb1->op == Add || mmb1->op == Sub || mmb1->op == Opp || mmb1->op == Div);
- g = g || (mmb2->op == Juxt || mmb2->op == Add || mmb2->op == Sub || mmb2->op == Opp || mmb2->op == Mult
- || mmb2->op == Div);
- if (f && g)
- sprintf(s, "(%s)/(%s)", s1, s2);
- else if (f)
- sprintf(s, "(%s)/%s", s1, s2);
- else if (g)
- sprintf(s, "%s/(%s)", s1, s2);
- else
- sprintf(s, "%s/%s", s1, s2);
- break;
- case Pow:
- f = (mmb1->op != Num && mmb1->op != Var);
- g = (mmb2->op != Num && mmb2->op != Var);
- if (f && g)
- sprintf(s, "(%s)^(%s)", s1, s2);
- else if (f)
- sprintf(s, "(%s)^%s", s1, s2);
- else if (g)
- sprintf(s, "%s^(%s)", s1, s2);
- else
- sprintf(s, "%s^%s", s1, s2);
- break;
- case Sqrt:
- g = (mmb2->op != Num && mmb2->op != Var && !g);
- if (g)
- sprintf(s, "sqrt(%s)", s2);
- else
- sprintf(s, "sqrt %s", s2);
- break;
- case NthRoot:
- f = (mmb1->op != Num && mmb1->op != Var);
- g = (mmb2->op != Num && mmb2->op != Var);
- if (f && g)
- sprintf(s, "(%s)#(%s)", s1, s2);
- else if (f)
- sprintf(s, "(%s)#%s", s1, s2);
- else if (g)
- sprintf(s, "%s#(%s)", s1, s2);
- else
- sprintf(s, "%s#%s", s1, s2);
- break;
- case E10:
- f = (mmb1->op != Num && mmb1->op != Var);
- g = (mmb2->op != Num && mmb2->op != Var);
- if (f && g)
- sprintf(s, "(%s)E(%s)", s1, s2);
- else if (f)
- sprintf(s, "(%s)E%s", s1, s2);
- else if (g)
- sprintf(s, "%sE(%s)", s1, s2);
- else
- sprintf(s, "%sE%s", s1, s2);
- break;
- case Ln:
- g = (mmb2->op != Num && mmb2->op != Var && !g);
- if (g)
- sprintf(s, "log(%s)", s2);
- else
- sprintf(s, "log %s", s2);
- break;
- case Exp:
- g = (mmb2->op != Num && mmb2->op != Var && !g);
- if (g)
- sprintf(s, "exp(%s)", s2);
- else
- sprintf(s, "exp %s", s2);
- break;
- case Sin:
- g = (mmb2->op != Num && mmb2->op != Var && !g);
- if (g)
- sprintf(s, "sin(%s)", s2);
- else
- sprintf(s, "sin %s", s2);
- break;
- case Cos:
- g = (mmb2->op != Num && mmb2->op != Var && !g);
- if (g)
- sprintf(s, "cos(%s)", s2);
- else
- sprintf(s, "cos %s", s2);
- break;
- case Tg:
- g = (mmb2->op != Num && mmb2->op != Var && !g);
- if (g)
- sprintf(s, "tan(%s)", s2);
- else
- sprintf(s, "tan %s", s2);
- break;
- case Atan:
- g = (mmb2->op != Num && mmb2->op != Var && !g);
- if (g)
- sprintf(s, "atan(%s)", s2);
- else
- sprintf(s, "atan %s", s2);
- break;
- case Asin:
- g = (mmb2->op != Num && mmb2->op != Var && !g);
- if (g)
- sprintf(s, "asin(%s)", s2);
- else
- sprintf(s, "asin %s", s2);
- break;
- case Acos:
- g = (mmb2->op != Num && mmb2->op != Var && !g);
- if (g)
- sprintf(s, "acos(%s)", s2);
- else
- sprintf(s, "acos %s", s2);
- break;
- case Abs:
- g = (mmb2->op != Num && mmb2->op != Var && !g);
- if (g)
- sprintf(s, "abs(%s)", s2);
- else
- sprintf(s, "abs %s", s2);
- break;
- case Fun:
- sprintf(s, "%s(%s)", pfunc->name, s2);
- break;
- default:
- return CopyStr("Error");
+ char*s=NULL,*s1=NULL,*s2=NULL;int n=10;signed char f=0,g=0;
+ if(op==Fun)if(strlen(pfunc->name)>4)n+=strlen(pfunc->name)-4;
+ if(mmb1!=NULL){s1=mmb1->Expr();n+=strlen(s1);f=IsFunction(mmb1->op);}
+ if(mmb2!=NULL){s2=mmb2->Expr();n+=strlen(s2);g=IsFunction(mmb2->op);}
+ s=new char[n];
+ switch(op){
+ case Num:return ValToStr(ValC);
+ case Var:return CopyStr(pvar->name);
+ case Juxt:sprintf(s,"%s , %s",s1,s2);break;
+ case Add:
+ f=f||(mmb1->op==Juxt);
+ g=g||(mmb2->op==Juxt);
+ if(f&&g)sprintf(s,"(%s)+(%s)",s1,s2);else
+ if(f)sprintf(s,"(%s)+%s",s1,s2);else
+ if(g)sprintf(s,"%s+(%s)",s1,s2);else
+ sprintf(s,"%s+%s",s1,s2);
+ break;
+ case Sub:
+ f=f||(mmb1->op==Juxt);
+ g=g||(mmb2->op==Juxt||mmb2->op==Add||mmb2->op==Sub);
+ if(f&&g)sprintf(s,"(%s)-(%s)",s1,s2);else
+ if(f)sprintf(s,"(%s)-%s",s1,s2);else
+ if(g)sprintf(s,"%s-(%s)",s1,s2);else
+ sprintf(s,"%s-%s",s1,s2);
+ break;
+ case Opp:
+ if(mmb2->op==Add||mmb2->op==Sub||mmb2->op==Juxt)sprintf(s,"-(%s)",s2);else
+ sprintf(s,"-%s",s2);
+ break;
+ case Mult:
+ f=f||(mmb1->op==Juxt||mmb1->op==Add||mmb1->op==Sub||mmb1->op==Opp||mmb1->op==Div);
+ g=g||(mmb2->op==Juxt||mmb2->op==Add||mmb2->op==Sub||mmb2->op==Opp);
+ if(f&&g)sprintf(s,"(%s)*(%s)",s1,s2);else
+ if(f)sprintf(s,"(%s)*%s",s1,s2);else
+ if(g)sprintf(s,"%s*(%s)",s1,s2);else
+ sprintf(s,"%s*%s",s1,s2);
+ break;
+ case Div:
+ f=f||(mmb1->op==Juxt||mmb1->op==Add||mmb1->op==Sub||mmb1->op==Opp||mmb1->op==Div);
+ g=g||(mmb2->op==Juxt||mmb2->op==Add||mmb2->op==Sub||mmb2->op==Opp||mmb2->op==Mult||mmb2->op==Div);
+ if(f&&g)sprintf(s,"(%s)/(%s)",s1,s2);else
+ if(f)sprintf(s,"(%s)/%s",s1,s2);else
+ if(g)sprintf(s,"%s/(%s)",s1,s2);else
+ sprintf(s,"%s/%s",s1,s2);
+ break;
+ case Pow:
+ f=(mmb1->op!=Num&&mmb1->op!=Var);
+ g=(mmb2->op!=Num&&mmb2->op!=Var);
+ if(f&&g)sprintf(s,"(%s)^(%s)",s1,s2);else
+ if(f)sprintf(s,"(%s)^%s",s1,s2);else
+ if(g)sprintf(s,"%s^(%s)",s1,s2);else
+ sprintf(s,"%s^%s",s1,s2);
+ break;
+ case Sqrt:
+ g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+ if(g)sprintf(s,"sqrt(%s)",s2);
+ else sprintf(s,"sqrt %s",s2);
+ break;
+ case NthRoot:
+ f=(mmb1->op!=Num&&mmb1->op!=Var);
+ g=(mmb2->op!=Num&&mmb2->op!=Var);
+ if(f&&g)sprintf(s,"(%s)#(%s)",s1,s2);else
+ if(f)sprintf(s,"(%s)#%s",s1,s2);else
+ if(g)sprintf(s,"%s#(%s)",s1,s2);else
+ sprintf(s,"%s#%s",s1,s2);
+ break;
+ case E10:
+ f=(mmb1->op!=Num&&mmb1->op!=Var);
+ g=(mmb2->op!=Num&&mmb2->op!=Var);
+ if(f&&g)sprintf(s,"(%s)E(%s)",s1,s2);else
+ if(f)sprintf(s,"(%s)E%s",s1,s2);else
+ if(g)sprintf(s,"%sE(%s)",s1,s2);else
+ sprintf(s,"%sE%s",s1,s2);
+ break;
+ case Ln:
+ g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+ if(g)sprintf(s,"log(%s)",s2);
+ else sprintf(s,"log %s",s2);
+ break;
+ case Exp:
+ g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+ if(g)sprintf(s,"exp(%s)",s2);
+ else sprintf(s,"exp %s",s2);
+ break;
+ case Sin:
+ g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+ if(g)sprintf(s,"sin(%s)",s2);
+ else sprintf(s,"sin %s",s2);
+ break;
+ case Cos:
+ g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+ if(g)sprintf(s,"cos(%s)",s2);
+ else sprintf(s,"cos %s",s2);
+ break;
+ case Tg:
+ g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+ if(g)sprintf(s,"tan(%s)",s2);
+ else sprintf(s,"tan %s",s2);
+ break;
+ case Atan:
+ g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+ if(g)sprintf(s,"atan(%s)",s2);
+ else sprintf(s,"atan %s",s2);
+ break;
+ case Asin:
+ g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+ if(g)sprintf(s,"asin(%s)",s2);
+ else sprintf(s,"asin %s",s2);
+ break;
+ case Acos:
+ g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+ if(g)sprintf(s,"acos(%s)",s2);
+ else sprintf(s,"acos %s",s2);
+ break;
+ case Abs:
+ g=(mmb2->op!=Num&&mmb2->op!=Var&&!g);
+ if(g)sprintf(s,"abs(%s)",s2);
+ else sprintf(s,"abs %s",s2);
+ break;
+ case Fun:
+ sprintf(s,"%s(%s)",pfunc->name,s2);
+ break;
+ default:return CopyStr("Error");
};
- if (s1 != NULL)
- delete[] s1;
- if (s2 != NULL)
- delete[] s2;
+ if(s1!=NULL)delete[] s1;if(s2!=NULL)delete[] s2;
return s;
}
-const double sqrtmaxfloat = sqrt(DBL_MAX);
-const double sqrtminfloat = sqrt(DBL_MIN);
-const double inveps = .1 / DBL_EPSILON;
-
-void Addition(double*& p)
-{
- if (*p == ErrVal || fabsl(*p) > sqrtmaxfloat) {
- *(--p) = ErrVal;
- return;
- };
- if (*(--p) == ErrVal || fabsl(*p) > sqrtmaxfloat) {
- *p = ErrVal;
- return;
- };
- *p += (*(p + 1));
-}
-void Soustraction(double*& p)
-{
- if (*p == ErrVal || fabsl(*p) > sqrtmaxfloat) {
- *(--p) = ErrVal;
- return;
- };
- if (*(--p) == ErrVal || fabsl(*p) > sqrtmaxfloat) {
- *p = ErrVal;
- return;
- };
- *p -= (*(p + 1));
-}
-void Multiplication(double*& p)
-{
- if (fabsl(*p) < sqrtminfloat) {
- *--p = 0;
- return;
- };
- if (*p == ErrVal || fabsl(*p) > sqrtmaxfloat) {
- *(--p) = ErrVal;
- return;
- };
- if (fabsl(*(--p)) < sqrtminfloat) {
- *p = 0;
- return;
- };
- if (*p == ErrVal || fabsl(*p) > sqrtmaxfloat) {
- *p = ErrVal;
- return;
- };
- *p *= (*(p + 1));
-}
-void Division(double*& p)
-{
- if (fabsl(*p) < sqrtminfloat || *p == ErrVal || fabsl(*p) > sqrtmaxfloat) {
- *(--p) = ErrVal;
- return;
- };
- if (fabsl(*(--p)) < sqrtminfloat)
- *p = 0;
- else if (*p == ErrVal || fabsl(*p) > sqrtmaxfloat) {
- *p = ErrVal;
- return;
- };
- *p /= (*(p + 1));
-}
-void Puissance(double*& p)
-{
- double v2 = *p--, v1 = *p;
- if (!v1) {
- *p = 0;
- return;
- };
- if (v2 == ErrVal || v1 == ErrVal || fabsl(v2 * logl(fabsl(v1))) > DBL_MAX_EXP) {
- *p = ErrVal;
- return;
- };
- *p = ((v1 > 0 || !fmodl(v2, 1)) ? powl(v1, v2) : ErrVal);
-}
-void RacineN(double*& p)
-{
- double v2 = *p--, v1 = *p;
- if (v1 == ErrVal || v2 == ErrVal || !v1 || v2 * logl(fabsl(v1)) < DBL_MIN_EXP) {
- *p = ErrVal;
- return;
- };
- if (v2 >= 0) {
- *p = powl(v2, 1 / v1);
- return;
- };
- *p = ((fabsl(fmodl(v1, 2)) == 1) ? -powl(-v2, 1 / v1) : ErrVal);
-}
-void Puiss10(double*& p)
-{
- if (fabsl(*p) < sqrtminfloat) {
- *(--p) = 0;
- return;
- };
- if (*p == ErrVal || fabsl(*p) > DBL_MAX_10_EXP) {
- *(--p) = ErrVal;
- return;
- };
- if (fabsl(*(--p)) < sqrtminfloat)
- *p = 0;
- else if (*p == ErrVal || fabsl(*p) > sqrtmaxfloat) {
- *p = ErrVal;
- return;
- };
- *p *= pow10l(*(p + 1));
-}
-void ArcTangente2(double*& p)
-{
- if (*p == ErrVal || fabsl(*p) > inveps) {
- *(--p) = ErrVal;
- return;
- };
- if (*(--p) == ErrVal || fabsl(*p) > inveps) {
- *p = ErrVal;
- return;
- };
- *p = (*p || *(p + 1) ? atan2(*p, *(p + 1)) : ErrVal);
-}
-void NextVal(double*&)
-{
-}
-void RFunc(double*&)
-{
-}
-void JuxtF(double*&)
-{
-}
-void Absolu(double*& p)
-{
- *p = ((*p == ErrVal) ? ErrVal : fabsl(*p));
-}
-void Oppose(double*& p)
-{
- *p = ((*p == ErrVal) ? ErrVal : -*p);
-}
-void ArcSinus(double*& p)
-{
- *p = ((*p == ErrVal || fabsl(*p) > 1) ? ErrVal : asinl(*p));
-}
-void ArcCosinus(double*& p)
-{
- *p = ((*p == ErrVal || fabsl(*p) > 1) ? ErrVal : acosl(*p));
-}
-void ArcTangente(double*& p)
-{
- *p = ((*p == ErrVal) ? ErrVal : atanl(*p));
-}
-void Logarithme(double*& p)
-{
- *p = ((*p == ErrVal || *p <= 0) ? ErrVal : logl(*p));
-}
-void Exponentielle(double*& p)
-{
- *p = ((*p == ErrVal || *p > DBL_MAX_EXP) ? ErrVal : expl(*p));
-}
-void Sinus(double*& p)
-{
- *p = ((*p == ErrVal || fabsl(*p) > inveps) ? ErrVal : sinl(*p));
-}
-void Tangente(double*& p)
-{
- *p = ((*p == ErrVal || fabsl(*p) > inveps) ? ErrVal : tanl(*p));
-}
-void Cosinus(double*& p)
-{
- *p = ((*p == ErrVal || fabsl(*p) > inveps) ? ErrVal : cosl(*p));
-}
-void Racine(double*& p)
-{
- *p = ((*p == ErrVal || *p > sqrtmaxfloat || *p < 0) ? ErrVal : sqrtl(*p));
-}
-void FonctionError(double*& p)
-{
- *p = ErrVal;
-}
-inline void ApplyRFunc(PRFunction rf, double*& p)
-{
- p -= rf->nvars - 1;
- *p = rf->Val(p);
-}
+const double sqrtmaxfloat=sqrt(DBL_MAX);
+const double sqrtminfloat=sqrt(DBL_MIN);
+const double inveps=.1/DBL_EPSILON;
+
+void Addition(double*&p)
+{if(*p==ErrVal||fabsl(*p)>sqrtmaxfloat){*(--p)=ErrVal;return;};
+if(*(--p)==ErrVal||fabsl(*p)>sqrtmaxfloat){*p=ErrVal;return;};
+*p+=(*(p+1));}
+void Soustraction(double*&p)
+{if(*p==ErrVal||fabsl(*p)>sqrtmaxfloat){*(--p)=ErrVal;return;};
+if(*(--p)==ErrVal||fabsl(*p)>sqrtmaxfloat){*p=ErrVal;return;};
+*p-=(*(p+1));}
+void Multiplication(double*&p)
+{if(fabsl(*p)<sqrtminfloat){*--p=0;return;};
+if(*p==ErrVal||fabsl(*p)>sqrtmaxfloat){*(--p)=ErrVal;return;};
+if(fabsl(*(--p))<sqrtminfloat){*p=0;return;};
+if(*p==ErrVal||fabsl(*p)>sqrtmaxfloat){*p=ErrVal;return;};
+*p*=(*(p+1));}
+void Division(double*&p)
+{if(fabsl(*p)<sqrtminfloat||*p==ErrVal||fabsl(*p)>sqrtmaxfloat)
+ {*(--p)=ErrVal;return;};
+if(fabsl(*(--p))<sqrtminfloat)*p=0;else if(*p==ErrVal||fabsl(*p)>sqrtmaxfloat)
+ {*p=ErrVal;return;};
+*p/=(*(p+1));}
+void Puissance(double*&p)
+{double v2=*p--,v1=*p;
+if(!v1){*p=0;return;};
+if(v2==ErrVal||v1==ErrVal||fabsl(v2*logl(fabsl(v1)))>DBL_MAX_EXP){*p=ErrVal;return;};
+*p=((v1>0||!fmodl(v2,1))?powl(v1,v2):ErrVal);}
+void RacineN(double*&p)
+{double v2=*p--,v1=*p;
+if(v1==ErrVal||v2==ErrVal||!v1||v2*logl(fabsl(v1))<DBL_MIN_EXP){*p=ErrVal;return;};
+if(v2>=0){*p=powl(v2,1/v1);return;};
+*p=((fabsl(fmodl(v1,2))==1)?-powl(-v2,1/v1):ErrVal);}
+void Puiss10(double*&p)
+{if(fabsl(*p)<sqrtminfloat){*(--p)=0;return;};
+if(*p==ErrVal||fabsl(*p)>DBL_MAX_10_EXP){*(--p)=ErrVal;return;};
+if(fabsl(*(--p))<sqrtminfloat)*p=0;else if(*p==ErrVal||fabsl(*p)>sqrtmaxfloat)
+{*p=ErrVal;return;};
+*p*=pow10l(*(p+1));}
+void ArcTangente2(double*&p)
+{if(*p==ErrVal||fabsl(*p)>inveps){*(--p)=ErrVal;return;};
+if(*(--p)==ErrVal||fabsl(*p)>inveps){*p=ErrVal;return;};
+*p=(*p||*(p+1)?atan2(*p,*(p+1)):ErrVal);}
+void NextVal(double*&){}
+void RFunc(double*&){}
+void JuxtF(double*&){}
+void Absolu(double*&p){*p=((*p==ErrVal)?ErrVal:fabsl(*p));}
+void Oppose(double*&p){*p=((*p==ErrVal)?ErrVal:-*p);}
+void ArcSinus(double*&p)
+{*p=((*p==ErrVal||fabsl(*p)>1)?ErrVal:asinl(*p));}
+void ArcCosinus(double*&p)
+{*p=((*p==ErrVal||fabsl(*p)>1)?ErrVal:acosl(*p));}
+void ArcTangente(double*&p)
+{*p=((*p==ErrVal)?ErrVal:atanl(*p));}
+void Logarithme(double*&p)
+{*p=((*p==ErrVal||*p<=0)?ErrVal:logl(*p));}
+void Exponentielle(double*&p)
+{*p=((*p==ErrVal||*p>DBL_MAX_EXP)?ErrVal:expl(*p));}
+void Sinus(double*&p)
+{*p=((*p==ErrVal||fabsl(*p)>inveps)?ErrVal:sinl(*p));}
+void Tangente(double*&p)
+{*p=((*p==ErrVal||fabsl(*p)>inveps)?ErrVal:tanl(*p));}
+void Cosinus(double*&p)
+{*p=((*p==ErrVal||fabsl(*p)>inveps)?ErrVal:cosl(*p));}
+void Racine(double*&p)
+{*p=((*p==ErrVal||*p>sqrtmaxfloat||*p<0)?ErrVal:sqrtl(*p));}
+void FonctionError(double*&p){*p=ErrVal;}
+inline void ApplyRFunc(PRFunction rf,double*&p)
+{p-=rf->nvars-1;*p=rf->Val(p);}
double ROperation::Val() const
{
- pfoncld* p1 = pinstr;
- double **p2 = pvals, *p3 = ppile - 1;
- PRFunction* p4 = pfuncpile;
- for (; *p1 != NULL; p1++)
- if (*p1 == &NextVal)
- *(++p3) = **(p2++);
- else if (*p1 == &RFunc)
- ApplyRFunc(*(p4++), p3);
- else
- (**p1)(p3);
+ pfoncld*p1=pinstr;double**p2=pvals,*p3=ppile-1;PRFunction*p4=pfuncpile;
+ for(;*p1!=NULL;p1++)
+ if(*p1==&NextVal)*(++p3)=**(p2++);else
+ if(*p1==&RFunc) ApplyRFunc(*(p4++),p3);
+ else (**p1)(p3);
return *p3;
}
-void BCDouble(pfoncld*& pf, pfoncld* pf1, pfoncld* pf2, double**& pv, double** pv1, double** pv2, double*& pp,
- double* pp1, double* pp2, RFunction**& prf, RFunction** prf1, RFunction** prf2, pfoncld f)
-{
- pfoncld *pf3, *pf4 = pf1;
- long n1, n2;
- for (n1 = 0; *pf4 != NULL; pf4++, n1++)
- ;
- for (n2 = 0, pf4 = pf2; *pf4 != NULL; pf4++, n2++)
- ;
- pf = new pfoncld[n1 + n2 + 2];
- for (pf3 = pf, pf4 = pf1; *pf4 != NULL; pf3++, pf4++)
- *pf3 = *pf4;
- for (pf4 = pf2; *pf4 != NULL; pf3++, pf4++)
- *pf3 = *pf4;
- *pf3++ = f;
- *pf3 = NULL; // delete[]pf1,pf2;
- double **pv3, **pv4 = pv1;
- for (n1 = 0; *pv4 != NULL; pv4++, n1++)
- ;
- for (n2 = 0, pv4 = pv2; *pv4 != NULL; pv4++, n2++)
- ;
- pv = new double*[n1 + n2 + 1];
- for (pv3 = pv, pv4 = pv1; *pv4 != NULL; pv3++, pv4++)
- *pv3 = *pv4;
- for (pv4 = pv2; *pv4 != NULL; pv3++, pv4++)
- *pv3 = *pv4;
- *pv3 = NULL; // delete[]pv1,pv2;
- double *pp3, *pp4 = pp1;
- for (n1 = 0; *pp4 != ErrVal; pp4++, n1++)
- ;
- for (n2 = 0, pp4 = pp2; *pp4 != ErrVal; pp4++, n2++)
- ;
- pp = new double[n1 + n2 + 1]; // Really need to add and not to take max(n1,n2) in case of Juxt operator
- for (pp3 = pp, pp4 = pp1; *pp4 != ErrVal; pp3++, pp4++)
- *pp3 = 0;
- for (pp4 = pp2; *pp4 != ErrVal; pp3++, pp4++)
- *pp3 = 0;
- *pp3 = ErrVal; // delete[]pp1,pp2;
- PRFunction *prf3, *prf4 = prf1;
- for (n1 = 0; *prf4 != NULL; prf4++, n1++)
- ;
- for (n2 = 0, prf4 = prf2; *prf4 != NULL; prf4++, n2++)
- ;
- prf = new PRFunction[n1 + n2 + 1];
- for (prf3 = prf, prf4 = prf1; *prf4 != NULL; prf3++, prf4++)
- *prf3 = *prf4;
- for (prf4 = prf2; *prf4 != NULL; prf3++, prf4++)
- *prf3 = *prf4;
- *prf3 = NULL; // delete[]prf1,prf2;
-}
-
-void BCSimple(pfoncld*& pf, pfoncld* pf1, double**& pv, double** pv1, double*& pp, double* pp1, RFunction**& prf,
- RFunction** prf1, pfoncld f)
+void BCDouble(pfoncld*&pf,pfoncld*pf1,pfoncld*pf2,
+ double**&pv,double**pv1,double**pv2,
+ double*&pp,double*pp1,double*pp2,
+ RFunction**&prf,RFunction**prf1,RFunction**prf2,
+ pfoncld f)
{
- pfoncld *pf3, *pf4 = pf1;
- long n;
- for (n = 0; *pf4 != NULL; pf4++, n++)
- ;
- pf = new pfoncld[n + 2];
- for (pf4 = pf1, pf3 = pf; *pf4 != NULL; pf3++, pf4++)
- *pf3 = *pf4;
- *pf3++ = f;
- *pf3 = NULL; // delete[]pf1;
- double **pv3, **pv4 = pv1;
- for (n = 0; *pv4 != NULL; pv4++, n++)
- ;
- pv = new double*[n + 1];
- for (pv3 = pv, pv4 = pv1; *pv4 != NULL; pv3++, pv4++)
- *pv3 = *pv4;
- *pv3 = NULL; // delete[]pv1;
- double *pp3, *pp4 = pp1;
- for (n = 0; *pp4 != ErrVal; pp4++, n++)
- ;
- pp = new double[n + 1];
- for (pp3 = pp, pp4 = pp1; *pp4 != ErrVal; pp3++, pp4++)
- *pp3 = 0;
- *pp3 = ErrVal; // delete[]pp1;
- RFunction **prf3, **prf4 = prf1;
- for (n = 0; *prf4 != NULL; prf4++, n++)
- ;
- prf = new RFunction*[n + 1];
- for (prf3 = prf, prf4 = prf1; *prf4 != NULL; prf3++, prf4++)
- *prf3 = *prf4;
- *prf3 = NULL; // delete[]prf1;
-}
-
-void BCFun(pfoncld*& pf, pfoncld* pf1, double**& pv, double** pv1, double*& pp, double* pp1, RFunction**& prf,
- RFunction** prf1, PRFunction rf)
-{
- pfoncld *pf3, *pf4 = pf1;
- long n;
- for (n = 0; *pf4 != NULL; pf4++, n++)
- ;
- pf = new pfoncld[n + 2];
- for (pf4 = pf1, pf3 = pf; *pf4 != NULL; pf3++, pf4++)
- *pf3 = *pf4;
- *pf3++ = &RFunc;
- *pf3 = NULL; // delete[]pf1;
- double **pv3, **pv4 = pv1;
- for (n = 0; *pv4 != NULL; pv4++, n++)
- ;
- pv = new double*[n + 1];
- for (pv3 = pv, pv4 = pv1; *pv4 != NULL; pv3++, pv4++)
- *pv3 = *pv4;
- *pv3 = NULL; // delete[]pv1;
- double *pp3, *pp4 = pp1;
- for (n = 0; *pp4 != ErrVal; pp4++, n++)
- ;
- pp = new double[n + 1];
- for (pp3 = pp, pp4 = pp1; *pp4 != ErrVal; pp3++, pp4++)
- *pp3 = 0;
- *pp3 = ErrVal; // delete[]pp1;
- PRFunction *prf3, *prf4 = prf1;
- for (n = 0; *prf4 != NULL; prf4++, n++)
- ;
- prf = new PRFunction[n + 2];
- for (prf4 = prf1, prf3 = prf; *prf4 != NULL; prf3++, prf4++)
- *prf3 = *prf4;
- *prf3++ = rf;
- *prf3 = NULL; // delete[]pf1;
+ pfoncld*pf3,*pf4=pf1;long n1,n2;
+ for(n1=0;*pf4!=NULL;pf4++,n1++)
+ ;
+ for(n2=0,pf4=pf2;*pf4!=NULL;pf4++,n2++)
+ ;
+ pf=new pfoncld[n1+n2+2];
+ for(pf3=pf,pf4=pf1;*pf4!=NULL;pf3++,pf4++)*pf3=*pf4;
+ for(pf4=pf2;*pf4!=NULL;pf3++,pf4++)*pf3=*pf4;
+ *pf3++=f;*pf3=NULL;//delete[]pf1,pf2;
+ double**pv3,**pv4=pv1;
+ for(n1=0;*pv4!=NULL;pv4++,n1++)
+ ;
+ for(n2=0,pv4=pv2;*pv4!=NULL;pv4++,n2++)
+ ;
+ pv=new double*[n1+n2+1];
+ for(pv3=pv,pv4=pv1;*pv4!=NULL;pv3++,pv4++)*pv3=*pv4;
+ for(pv4=pv2;*pv4!=NULL;pv3++,pv4++)*pv3=*pv4;
+ *pv3=NULL;//delete[]pv1,pv2;
+ double*pp3,*pp4=pp1;
+ for(n1=0;*pp4!=ErrVal;pp4++,n1++)
+ ;
+ for(n2=0,pp4=pp2;*pp4!=ErrVal;pp4++,n2++)
+ ;
+ pp=new double[n1+n2+1]; // Really need to add and not to take max(n1,n2) in case of Juxt operator
+ for(pp3=pp,pp4=pp1;*pp4!=ErrVal;pp3++,pp4++)*pp3=0;
+ for(pp4=pp2;*pp4!=ErrVal;pp3++,pp4++)*pp3=0;
+ *pp3=ErrVal;//delete[]pp1,pp2;
+ PRFunction*prf3,*prf4=prf1;
+ for(n1=0;*prf4!=NULL;prf4++,n1++)
+ ;
+ for(n2=0,prf4=prf2;*prf4!=NULL;prf4++,n2++)
+ ;
+ prf=new PRFunction[n1+n2+1];
+ for(prf3=prf,prf4=prf1;*prf4!=NULL;prf3++,prf4++)*prf3=*prf4;
+ for(prf4=prf2;*prf4!=NULL;prf3++,prf4++)*prf3=*prf4;
+ *prf3=NULL;//delete[]prf1,prf2;
+}
+
+void BCSimple(pfoncld*&pf,pfoncld*pf1,double**&pv,double**pv1,
+ double*&pp,double*pp1,RFunction**&prf,RFunction**prf1,pfoncld f)
+{
+ pfoncld*pf3,*pf4=pf1;long n;
+ for(n=0;*pf4!=NULL;pf4++,n++);
+ pf=new pfoncld[n+2];
+ for(pf4=pf1,pf3=pf;*pf4!=NULL;pf3++,pf4++)*pf3=*pf4;
+ *pf3++=f;*pf3=NULL;//delete[]pf1;
+ double**pv3,**pv4=pv1;
+ for(n=0;*pv4!=NULL;pv4++,n++);
+ pv=new double*[n+1];
+ for(pv3=pv,pv4=pv1;*pv4!=NULL;pv3++,pv4++)*pv3=*pv4;
+ *pv3=NULL;//delete[]pv1;
+ double*pp3,*pp4=pp1;
+ for(n=0;*pp4!=ErrVal;pp4++,n++);
+ pp=new double[n+1];
+ for(pp3=pp,pp4=pp1;*pp4!=ErrVal;pp3++,pp4++)*pp3=0;
+ *pp3=ErrVal;//delete[]pp1;
+ RFunction**prf3,**prf4=prf1;
+ for(n=0;*prf4!=NULL;prf4++,n++);
+ prf=new RFunction*[n+1];
+ for(prf3=prf,prf4=prf1;*prf4!=NULL;prf3++,prf4++)*prf3=*prf4;
+ *prf3=NULL;//delete[]prf1;
+}
+
+void BCFun(pfoncld*&pf,pfoncld*pf1,double**&pv,double**pv1,
+ double*&pp,double*pp1,RFunction**&prf,RFunction**prf1,PRFunction rf)
+{
+ pfoncld*pf3,*pf4=pf1;long n;
+ for(n=0;*pf4!=NULL;pf4++,n++);
+ pf=new pfoncld[n+2];
+ for(pf4=pf1,pf3=pf;*pf4!=NULL;pf3++,pf4++)*pf3=*pf4;
+ *pf3++=&RFunc;*pf3=NULL;//delete[]pf1;
+ double**pv3,**pv4=pv1;
+ for(n=0;*pv4!=NULL;pv4++,n++);
+ pv=new double*[n+1];
+ for(pv3=pv,pv4=pv1;*pv4!=NULL;pv3++,pv4++)*pv3=*pv4;
+ *pv3=NULL;//delete[]pv1;
+ double*pp3,*pp4=pp1;
+ for(n=0;*pp4!=ErrVal;pp4++,n++);
+ pp=new double[n+1];
+ for(pp3=pp,pp4=pp1;*pp4!=ErrVal;pp3++,pp4++)*pp3=0;
+ *pp3=ErrVal;//delete[]pp1;
+ PRFunction*prf3,*prf4=prf1;
+ for(n=0;*prf4!=NULL;prf4++,n++);
+ prf=new PRFunction[n+2];
+ for(prf4=prf1,prf3=prf;*prf4!=NULL;prf3++,prf4++)*prf3=*prf4;
+ *prf3++=rf;*prf3=NULL;//delete[]pf1;
}
void ROperation::BuildCode()
{
// if(mmb1!=NULL)mmb1->BuildCode();if(mmb2!=NULL)mmb2->BuildCode();
- if (pinstr != NULL) {
- delete[] pinstr;
- pinstr = NULL;
- }
- if (pvals != NULL) {
- delete[] pvals;
- pvals = NULL;
- } // does not delete pvals[0] in case it was to be deleted... (no way to know)
- if (ppile != NULL) {
- delete[] ppile;
- ppile = NULL;
- }
- if (pfuncpile != NULL) {
- delete[] pfuncpile;
- pfuncpile = NULL;
- }
- switch (op) {
- case ErrOp:
- pinstr = new pfoncld[2];
- pinstr[0] = &NextVal;
- pinstr[1] = NULL;
- pvals = new double*[2];
- pvals[0] = new double(ErrVal);
- pvals[1] = NULL;
- ppile = new double[2];
- ppile[0] = 0;
- ppile[1] = ErrVal;
- pfuncpile = new RFunction*[1];
- pfuncpile[0] = NULL;
- break;
- case Num:
- pinstr = new pfoncld[2];
- pinstr[0] = &NextVal;
- pinstr[1] = NULL;
- pvals = new double*[2];
- pvals[0] = new double(ValC);
- pvals[1] = NULL;
- ppile = new double[2];
- ppile[0] = 0;
- ppile[1] = ErrVal;
- pfuncpile = new RFunction*[1];
- pfuncpile[0] = NULL;
- break;
- case Var:
- pinstr = new pfoncld[2];
- pinstr[0] = &NextVal;
- pinstr[1] = NULL;
- pvals = new double*[2];
- pvals[0] = pvarval;
- pvals[1] = NULL;
- ppile = new double[2];
- ppile[0] = 0;
- ppile[1] = ErrVal;
- pfuncpile = new RFunction*[1];
- pfuncpile[0] = NULL;
- break;
- case Juxt:
- BCDouble(pinstr,
- mmb1->pinstr,
- mmb2->pinstr,
- pvals,
- mmb1->pvals,
- mmb2->pvals,
- ppile,
- mmb1->ppile,
- mmb2->ppile,
- pfuncpile,
- mmb1->pfuncpile,
- mmb2->pfuncpile,
- &JuxtF);
- break;
- case Add:
- BCDouble(pinstr,
- mmb1->pinstr,
- mmb2->pinstr,
- pvals,
- mmb1->pvals,
- mmb2->pvals,
- ppile,
- mmb1->ppile,
- mmb2->ppile,
- pfuncpile,
- mmb1->pfuncpile,
- mmb2->pfuncpile,
- &Addition);
- break;
- case Sub:
- BCDouble(pinstr,
- mmb1->pinstr,
- mmb2->pinstr,
- pvals,
- mmb1->pvals,
- mmb2->pvals,
- ppile,
- mmb1->ppile,
- mmb2->ppile,
- pfuncpile,
- mmb1->pfuncpile,
- mmb2->pfuncpile,
- &Soustraction);
- break;
- case Mult:
- BCDouble(pinstr,
- mmb1->pinstr,
- mmb2->pinstr,
- pvals,
- mmb1->pvals,
- mmb2->pvals,
- ppile,
- mmb1->ppile,
- mmb2->ppile,
- pfuncpile,
- mmb1->pfuncpile,
- mmb2->pfuncpile,
- &Multiplication);
- break;
- case Div:
- BCDouble(pinstr,
- mmb1->pinstr,
- mmb2->pinstr,
- pvals,
- mmb1->pvals,
- mmb2->pvals,
- ppile,
- mmb1->ppile,
- mmb2->ppile,
- pfuncpile,
- mmb1->pfuncpile,
- mmb2->pfuncpile,
- &Division);
- break;
- case Pow:
- BCDouble(pinstr,
- mmb1->pinstr,
- mmb2->pinstr,
- pvals,
- mmb1->pvals,
- mmb2->pvals,
- ppile,
- mmb1->ppile,
- mmb2->ppile,
- pfuncpile,
- mmb1->pfuncpile,
- mmb2->pfuncpile,
- &Puissance);
- break;
- case NthRoot:
- BCDouble(pinstr,
- mmb1->pinstr,
- mmb2->pinstr,
- pvals,
- mmb1->pvals,
- mmb2->pvals,
- ppile,
- mmb1->ppile,
- mmb2->ppile,
- pfuncpile,
- mmb1->pfuncpile,
- mmb2->pfuncpile,
- &RacineN);
- break;
- case E10:
- BCDouble(pinstr,
- mmb1->pinstr,
- mmb2->pinstr,
- pvals,
- mmb1->pvals,
- mmb2->pvals,
- ppile,
- mmb1->ppile,
- mmb2->ppile,
- pfuncpile,
- mmb1->pfuncpile,
- mmb2->pfuncpile,
- &Puiss10);
- break;
- case Opp:
- BCSimple(pinstr, mmb2->pinstr, pvals, mmb2->pvals, ppile, mmb2->ppile, pfuncpile, mmb2->pfuncpile, &Oppose);
- break;
- case Sin:
- BCSimple(pinstr, mmb2->pinstr, pvals, mmb2->pvals, ppile, mmb2->ppile, pfuncpile, mmb2->pfuncpile, &Sinus);
- break;
- case Sqrt:
- BCSimple(pinstr, mmb2->pinstr, pvals, mmb2->pvals, ppile, mmb2->ppile, pfuncpile, mmb2->pfuncpile, &Racine);
- break;
- case Ln:
- BCSimple(pinstr, mmb2->pinstr, pvals, mmb2->pvals, ppile, mmb2->ppile, pfuncpile, mmb2->pfuncpile, &Logarithme);
- break;
- case Exp:
- BCSimple(
- pinstr, mmb2->pinstr, pvals, mmb2->pvals, ppile, mmb2->ppile, pfuncpile, mmb2->pfuncpile, &Exponentielle);
- break;
- case Cos:
- BCSimple(pinstr, mmb2->pinstr, pvals, mmb2->pvals, ppile, mmb2->ppile, pfuncpile, mmb2->pfuncpile, &Cosinus);
- break;
- case Tg:
- BCSimple(pinstr, mmb2->pinstr, pvals, mmb2->pvals, ppile, mmb2->ppile, pfuncpile, mmb2->pfuncpile, &Tangente);
- break;
- case Atan:
- BCSimple(pinstr,
- mmb2->pinstr,
- pvals,
- mmb2->pvals,
- ppile,
- mmb2->ppile,
- pfuncpile,
- mmb2->pfuncpile,
- (mmb2->NMembers() > 1 ? &ArcTangente2 : &ArcTangente));
- break;
- case Asin:
- BCSimple(pinstr, mmb2->pinstr, pvals, mmb2->pvals, ppile, mmb2->ppile, pfuncpile, mmb2->pfuncpile, &ArcSinus);
- break;
- case Acos:
- BCSimple(pinstr, mmb2->pinstr, pvals, mmb2->pvals, ppile, mmb2->ppile, pfuncpile, mmb2->pfuncpile, &ArcCosinus);
- break;
- case Abs:
- BCSimple(pinstr, mmb2->pinstr, pvals, mmb2->pvals, ppile, mmb2->ppile, pfuncpile, mmb2->pfuncpile, &Absolu);
- break;
- case Fun:
- BCFun(pinstr, mmb2->pinstr, pvals, mmb2->pvals, ppile, mmb2->ppile, pfuncpile, mmb2->pfuncpile, pfunc);
- break;
- default:
- BCSimple(
- pinstr, mmb2->pinstr, pvals, mmb2->pvals, ppile, mmb2->ppile, pfuncpile, mmb2->pfuncpile, &FonctionError);
+ if(pinstr!=NULL){delete[]pinstr;pinstr=NULL;}
+ if(pvals!=NULL){delete[]pvals;pvals=NULL;}//does not delete pvals[0] in case it was to be deleted... (no way to know)
+ if(ppile!=NULL){delete[]ppile;ppile=NULL;}
+ if(pfuncpile!=NULL){delete[]pfuncpile;pfuncpile=NULL;}
+ switch(op){
+ case ErrOp:pinstr=new pfoncld[2];pinstr[0]=&NextVal;pinstr[1]=NULL;
+ pvals=new double*[2];pvals[0]=new double(ErrVal);pvals[1]=NULL;
+ ppile=new double[2];ppile[0]=0;ppile[1]=ErrVal;
+ pfuncpile=new RFunction*[1];pfuncpile[0]=NULL;
+ break;
+ case Num:pinstr=new pfoncld[2];pinstr[0]=&NextVal;pinstr[1]=NULL;
+ pvals=new double*[2];pvals[0]=new double(ValC);pvals[1]=NULL;
+ ppile=new double[2];ppile[0]=0;ppile[1]=ErrVal;
+ pfuncpile=new RFunction*[1];pfuncpile[0]=NULL;break;
+ case Var:pinstr=new pfoncld[2];pinstr[0]=&NextVal;pinstr[1]=NULL;
+ pvals=new double*[2];pvals[0]=pvarval;pvals[1]=NULL;
+ ppile=new double[2];ppile[0]=0;ppile[1]=ErrVal;
+ pfuncpile=new RFunction*[1];pfuncpile[0]=NULL;break;
+ case Juxt:BCDouble(pinstr,mmb1->pinstr,mmb2->pinstr,
+ pvals,mmb1->pvals,mmb2->pvals,ppile,mmb1->ppile,mmb2->ppile,pfuncpile,mmb1->pfuncpile,mmb2->pfuncpile,&JuxtF);
+ break;case Add:BCDouble(pinstr,mmb1->pinstr,mmb2->pinstr,
+ pvals,mmb1->pvals,mmb2->pvals,ppile,mmb1->ppile,mmb2->ppile,pfuncpile,mmb1->pfuncpile,mmb2->pfuncpile,&Addition);
+ break;case Sub:BCDouble(pinstr,mmb1->pinstr,mmb2->pinstr,
+ pvals,mmb1->pvals,mmb2->pvals,ppile,mmb1->ppile,mmb2->ppile,pfuncpile,mmb1->pfuncpile,mmb2->pfuncpile,&Soustraction);
+ break;case Mult:BCDouble(pinstr,mmb1->pinstr,mmb2->pinstr,
+ pvals,mmb1->pvals,mmb2->pvals,ppile,mmb1->ppile,mmb2->ppile,pfuncpile,mmb1->pfuncpile,mmb2->pfuncpile,&Multiplication);
+ break;case Div:BCDouble(pinstr,mmb1->pinstr,mmb2->pinstr,
+ pvals,mmb1->pvals,mmb2->pvals,ppile,mmb1->ppile,mmb2->ppile,pfuncpile,mmb1->pfuncpile,mmb2->pfuncpile,&Division);
+ break;case Pow:BCDouble(pinstr,mmb1->pinstr,mmb2->pinstr,
+ pvals,mmb1->pvals,mmb2->pvals,ppile,mmb1->ppile,mmb2->ppile,pfuncpile,mmb1->pfuncpile,mmb2->pfuncpile,&Puissance);
+ break;case NthRoot:BCDouble(pinstr,mmb1->pinstr,mmb2->pinstr,
+ pvals,mmb1->pvals,mmb2->pvals,ppile,mmb1->ppile,mmb2->ppile,pfuncpile,mmb1->pfuncpile,mmb2->pfuncpile,&RacineN);
+ break;case E10:BCDouble(pinstr,mmb1->pinstr,mmb2->pinstr,
+ pvals,mmb1->pvals,mmb2->pvals,ppile,mmb1->ppile,mmb2->ppile,pfuncpile,mmb1->pfuncpile,mmb2->pfuncpile,&Puiss10);
+ break;case Opp:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+ ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Oppose);
+ break;case Sin:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+ ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Sinus);
+ break;case Sqrt:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+ ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Racine);
+ break;case Ln:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+ ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Logarithme);
+ break;case Exp:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+ ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Exponentielle);
+ break;case Cos:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+ ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Cosinus);
+ break;case Tg:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+ ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Tangente);
+ break;case Atan:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+ ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,(mmb2->NMembers()>1?&ArcTangente2:&ArcTangente));
+ break;case Asin:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+ ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&ArcSinus);
+ break;case Acos:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+ ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&ArcCosinus);
+ break;case Abs:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+ ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&Absolu);
+ break;case Fun:BCFun(pinstr,mmb2->pinstr,pvals,mmb2->pvals,ppile,
+ mmb2->ppile,pfuncpile,mmb2->pfuncpile,pfunc);
+ break;default:BCSimple(pinstr,mmb2->pinstr,pvals,mmb2->pvals,
+ ppile,mmb2->ppile,pfuncpile,mmb2->pfuncpile,&FonctionError);
}
}
+
+// clang-format on
diff --git a/dune/stuff/functions/expression/mathexpr.hh b/dune/stuff/functions/expression/mathexpr.hh
index c01e038c8516fff72efa4fd2256fd8a2a9d2e1a2..05166ae5d21ee585f4813c6303b962b8c7d1c6fc 100644
--- a/dune/stuff/functions/expression/mathexpr.hh
+++ b/dune/stuff/functions/expression/mathexpr.hh
@@ -19,14 +19,16 @@ This software comes with absolutely no warranty.
*/
+// clang-format off
+
#ifndef DUNE_STUFF_FUNCTION_NONPARAMETRIC_EXPRESSION_MATHEXPRESSION_HH
#define DUNE_STUFF_FUNCTION_NONPARAMETRIC_EXPRESSION_MATHEXPRESSION_HH
-#include <string.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-#include <float.h>
+#include<string.h>
+#include<stdio.h>
+#include<stdlib.h>
+#include<math.h>
+#include<float.h>
// Compatibility with long double-typed functions
#define atanl atan
@@ -35,7 +37,7 @@ This software comes with absolutely no warranty.
#define expl exp
#define logl log
#define powl pow
-#define pow10l(x) pow(10, x)
+#define pow10l(x) pow(10,x)
#define fabsl fabs
#define cosl cos
#define sinl sin
@@ -45,12 +47,11 @@ This software comes with absolutely no warranty.
// Warning : if ints are short, everything will fail with strings longer than 32767 chars
-const double ErrVal = DBL_MAX;
+const double ErrVal=DBL_MAX;
-// Class definitions for operations
+//Class definitions for operations
-class RVar
-{
+class RVar{
public:
RVar()
{
@@ -66,88 +67,54 @@ public:
friend int operator==(const RVar&, const RVar&);
- char* name;
- double* pval;
+ char*name;double*pval;
};
typedef RVar* PRVar;
-enum ROperator
-{
- ErrOp,
- Juxt,
- Num,
- Var,
- Add,
- Sub,
- Opp,
- Mult,
- Div,
- Pow,
- Sqrt,
- NthRoot,
- Abs,
- Sin,
- Cos,
- Tg,
- Ln,
- Exp,
- Acos,
- Asin,
- Atan,
- E10,
- Fun
-};
+enum ROperator{ErrOp,Juxt,Num,Var,Add,Sub,Opp,Mult,Div,Pow,Sqrt,
+ NthRoot,Abs,Sin,Cos,Tg,Ln,Exp,Acos,Asin,Atan,E10,Fun};
-typedef void((*pfoncld)(double*&));
+typedef void ((*pfoncld)(double*&));
class ROperation;
typedef ROperation* PROperation;
class RFunction;
typedef RFunction* PRFunction;
-class ROperation
-{
- pfoncld* pinstr;
- double** pvals;
- double* ppile;
- RFunction** pfuncpile;
+class ROperation{
+ pfoncld*pinstr;double**pvals;double*ppile;RFunction**pfuncpile;
mutable signed char containfuncflag;
void BuildCode();
void Destroy();
-
-public:
+ public:
ROperator op;
- PROperation mmb1, mmb2;
- double ValC;
- const RVar* pvar;
- double* pvarval;
+ PROperation mmb1,mmb2;
+ double ValC;const RVar* pvar;double*pvarval;
RFunction* pfunc;
ROperation();
ROperation(const ROperation&);
ROperation(double);
ROperation(const RVar&);
- ROperation(const char* sp, int nvarp = 0, PRVar* ppvarp = NULL, int nfuncp = 0, PRFunction* ppfuncp = NULL);
+ ROperation(const char*sp,int nvarp=0,PRVar*ppvarp=NULL,int nfuncp=0,PRFunction*ppfuncp=NULL);
~ROperation();
double Val() const;
signed char ContainVar(const RVar&) const;
signed char ContainFunc(const RFunction&) const;
signed char ContainFuncNoRec(const RFunction&) const; // No recursive test on subfunctions
- ROperation NthMember(int) const;
- int NMembers() const;
- signed char HasError(const ROperation* = NULL) const;
+ ROperation NthMember(int) const;int NMembers() const;
+ signed char HasError(const ROperation* =NULL) const;
ROperation& operator=(const ROperation&);
- friend int operator==(const ROperation&, const double);
- friend int operator==(const ROperation&, const ROperation&);
- friend int operator!=(const ROperation&, const ROperation&);
- ROperation operator+() const;
- ROperation operator-() const;
- friend ROperation operator, (const ROperation&, const ROperation&);
- friend ROperation operator+(const ROperation&, const ROperation&);
- friend ROperation operator-(const ROperation&, const ROperation&);
- friend ROperation operator*(const ROperation&, const ROperation&);
- friend ROperation operator/(const ROperation&, const ROperation&);
- friend ROperation operator^(const ROperation&, const ROperation&); // Caution: wrong associativity and precedence
+ friend int operator==(const ROperation& ,const double);
+ friend int operator==(const ROperation& ,const ROperation&);
+ friend int operator!=(const ROperation& ,const ROperation&);
+ ROperation operator+() const;ROperation operator-() const;
+ friend ROperation operator,(const ROperation&,const ROperation&);
+ friend ROperation operator+(const ROperation&,const ROperation&);
+ friend ROperation operator-(const ROperation&,const ROperation&);
+ friend ROperation operator*(const ROperation&,const ROperation&);
+ friend ROperation operator/(const ROperation&,const ROperation&);
+ friend ROperation operator^(const ROperation&,const ROperation&); // Caution: wrong associativity and precedence
friend ROperation sqrt(const ROperation&);
friend ROperation abs(const ROperation&);
friend ROperation sin(const ROperation&);
@@ -158,42 +125,40 @@ public:
friend ROperation acos(const ROperation&);
friend ROperation asin(const ROperation&);
friend ROperation atan(const ROperation&);
- friend ROperation ApplyOperator(int, ROperation**, ROperation (*)(const ROperation&, const ROperation&));
+ friend ROperation ApplyOperator(int,ROperation**,ROperation (*)(const ROperation&,const ROperation&));
ROperation Diff(const RVar&) const; // Differentiate w.r.t a variable
char* Expr() const;
- ROperation Substitute(const RVar&, const ROperation&) const;
+ ROperation Substitute(const RVar&,const ROperation&) const;
};
-class RFunction
-{
- double* buf;
-
+class RFunction{
+ double*buf;
public:
signed char type;
- double((*pfuncval)(double));
- ROperation op;
- int nvars;
- RVar** ppvar;
- char* name;
+ double ((*pfuncval)(double));
+ ROperation op;int nvars;RVar** ppvar;
+ char*name;
RFunction();
- RFunction(double((*)(double)));
- RFunction(const ROperation& opp, RVar* pvarp);
- RFunction(const ROperation& opp, int nvarsp, RVar** ppvarp);
+ RFunction(double ((*)(double)));
+ RFunction(const ROperation& opp,RVar* pvarp);
+ RFunction(const ROperation& opp,int nvarsp,RVar**ppvarp);
RFunction(const RFunction&);
~RFunction();
RFunction& operator=(const RFunction&);
- void SetName(const char* s);
+ void SetName(const char*s);
double Val(double) const;
double Val(double*) const;
- friend int operator==(const RFunction&, const RFunction&);
+ friend int operator==(const RFunction&,const RFunction&);
ROperation operator()(const ROperation&);
};
-char* MidStr(const char* s, int i1, int i2);
-char* CopyStr(const char* s);
-char* InsStr(const char* s, int n, char c);
-signed char EqStr(const char* s, const char* s2);
-signed char CompStr(const char* s, int n, const char* s2);
-char* DelStr(const char* s, int n);
+char* MidStr(const char*s,int i1,int i2);
+char* CopyStr(const char*s);
+char* InsStr(const char*s,int n,char c);
+signed char EqStr(const char*s,const char*s2);
+signed char CompStr(const char*s,int n,const char*s2);
+char* DelStr(const char*s,int n);
#endif // DUNE_STUFF_FUNCTION_NONPARAMETRIC_EXPRESSION_MATHEXPRESSION_HH
+
+// clang-format on
diff --git a/dune/stuff/functions/indicator.hh b/dune/stuff/functions/indicator.hh
new file mode 100644
index 0000000000000000000000000000000000000000..5fe128ba4c6e611bf599ce314959bed10b0769e1
--- /dev/null
+++ b/dune/stuff/functions/indicator.hh
@@ -0,0 +1,193 @@
+// This file is part of the dune-stuff project:
+// https://github.com/wwu-numerik/dune-stuff
+// Copyright holders: Rene Milk, Felix Schindler
+// License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+
+#ifndef DUNE_STUFF_FUNCTIONS_INDICATOR_HH
+#define DUNE_STUFF_FUNCTIONS_INDICATOR_HH
+
+#include <memory>
+#include <vector>
+#include <utility>
+
+#include <dune/stuff/common/type_utils.hh>
+
+#include <dune/stuff/common/fvector.hh>
+#include <dune/stuff/common/string.hh>
+#include <dune/stuff/common/memory.hh>
+#include <dune/stuff/common/configuration.hh>
+
+#include "interfaces.hh"
+
+namespace Dune {
+namespace Stuff {
+namespace Functions {
+
+
+template <class E, class D, size_t d, class R, size_t r, size_t rC = 1>
+class Indicator : public LocalizableFunctionInterface<E, D, d, R, r, rC>
+{
+ Indicator()
+ {
+ static_assert(AlwaysFalse<E>::value, "Not available for these dimensions!");
+ }
+};
+
+
+template <class E, class D, size_t d, class R>
+class Indicator<E, D, d, R, 1> : public LocalizableFunctionInterface<E, D, d, R, 1>
+{
+ typedef LocalizableFunctionInterface<E, D, d, R, 1> BaseType;
+ typedef Indicator<E, D, d, R, 1> ThisType;
+
+ class Localfunction : public LocalfunctionInterface<E, D, d, R, 1>
+ {
+ typedef LocalfunctionInterface<E, D, d, R, 1> InterfaceType;
+
+ public:
+ using typename InterfaceType::EntityType;
+ using typename InterfaceType::DomainType;
+ using typename InterfaceType::RangeType;
+ using typename InterfaceType::JacobianRangeType;
+
+ Localfunction(const EntityType& entity, const RangeType& value)
+ : InterfaceType(entity)
+ , value_(value)
+ {
+ }
+
+ virtual size_t order() const override final
+ {
+ return 0;
+ }
+
+ virtual void evaluate(const DomainType& xx, RangeType& ret) const override final
+ {
+ assert(this->is_a_valid_point(xx));
+ ret = value_;
+ }
+
+ virtual void jacobian(const DomainType& xx, JacobianRangeType& ret) const override final
+ {
+ assert(this->is_a_valid_point(xx));
+ ret *= 0.0;
+ }
+
+ private:
+ const RangeType value_;
+ }; // class Localfunction
+
+public:
+ using typename BaseType::EntityType;
+ using typename BaseType::DomainFieldType;
+ using typename BaseType::DomainType;
+ using typename BaseType::RangeType;
+ using typename BaseType::RangeFieldType;
+ using typename BaseType::LocalfunctionType;
+
+ static const bool available = true;
+
+ static std::string static_id()
+ {
+ return BaseType::static_id() + ".indicator";
+ }
+
+ static Common::Configuration default_config(const std::string sub_name = "")
+ {
+ Common::Configuration cfg;
+ cfg["name"] = static_id();
+ if (d == 1)
+ cfg["0.domain"] = "[0.25 0.75]";
+ else if (d == 2)
+ cfg["0.domain"] = "[0.25 0.75; 0.25 0.75]";
+ else if (d == 3)
+ cfg["0.domain"] = "[0.25 0.75; 0.25 0.75; 0.25 0.75]";
+ else
+ DUNE_THROW(NotImplemented, "Indeed!");
+ cfg["0.value"] = "1";
+ if (sub_name.empty())
+ return cfg;
+ else {
+ Common::Configuration tmp;
+ tmp.add(cfg, sub_name);
+ return tmp;
+ }
+ } // ... default_config(...)
+
+ static std::unique_ptr<ThisType> create(const Common::Configuration config = default_config(),
+ const std::string sub_name = static_id())
+ {
+ const Common::Configuration cfg = config.has_sub(sub_name) ? config.sub(sub_name) : config;
+ const Common::Configuration def_cfg = default_config();
+ std::vector<std::tuple<DomainType, DomainType, RangeFieldType>> values;
+ DomainType tmp_lower;
+ DomainType tmp_upper;
+ size_t cc = 0;
+ while (cfg.has_sub(DSC::toString(cc))) {
+ const Stuff::Common::Configuration local_cfg = cfg.sub(DSC::toString(cc));
+ if (local_cfg.has_key("domain") && local_cfg.has_key("value")) {
+ auto domains = local_cfg.get<FieldMatrix<DomainFieldType, d, 2>>("domain");
+ for (size_t dd = 0; dd < d; ++dd) {
+ tmp_lower[dd] = domains[dd][0];
+ tmp_upper[dd] = domains[dd][1];
+ }
+ auto val = local_cfg.get<RangeFieldType>("value");
+ values.emplace_back(tmp_lower, tmp_upper, val);
+ } else
+ break;
+ ++cc;
+ }
+ return Common::make_unique<ThisType>(values, cfg.get("name", def_cfg.get<std::string>("name")));
+ } // ... create(...)
+
+ Indicator(const std::vector<std::tuple<DomainType, DomainType, R>>& values, const std::string name_in = "indicator")
+ : values_(values)
+ , name_(name_in)
+ {
+ }
+
+ Indicator(const std::vector<std::pair<std::pair<Common::FieldVector<D, d>, Common::FieldVector<D, d>>, R>>& values,
+ const std::string name_in = "indicator")
+ : values_(convert(values))
+ , name_(name_in)
+ {
+ }
+
+ virtual ~Indicator()
+ {
+ }
+
+ virtual std::string name() const override final
+ {
+ return name_;
+ }
+
+ virtual std::unique_ptr<LocalfunctionType> local_function(const EntityType& entity) const override final
+ {
+ const auto center = entity.geometry().center();
+ for (const auto& element : values_)
+ if (Common::FloatCmp::le(std::get<0>(element), center) && Common::FloatCmp::lt(center, std::get<1>(element)))
+ return Common::make_unique<Localfunction>(entity, std::get<2>(element));
+ return Common::make_unique<Localfunction>(entity, 0.0);
+ } // ... local_function(...)
+
+private:
+ static std::vector<std::tuple<DomainType, DomainType, R>>
+ convert(const std::vector<std::pair<std::pair<Common::FieldVector<D, d>, Common::FieldVector<D, d>>, R>>& values)
+ {
+ std::vector<std::tuple<DomainType, DomainType, R>> ret;
+ for (const auto& element : values)
+ ret.emplace_back(element.first.first, element.first.second, element.second);
+ return ret;
+ } // convert(...)
+
+ const std::vector<std::tuple<DomainType, DomainType, R>> values_;
+ const std::string name_;
+}; // class Indicator
+
+
+} // namespace Functions
+} // namespace Stuff
+} // namespace Dune
+
+#endif // DUNE_STUFF_FUNCTIONS_INDICATOR_HH
diff --git a/dune/stuff/functions/random_ellipsoids_function.hh b/dune/stuff/functions/random_ellipsoids_function.hh
index b9edc8ae1fd86519e9e5899ff0248bc9b50aba5b..a0036cb174b1836340f8f4c195a68def415ebffc 100644
--- a/dune/stuff/functions/random_ellipsoids_function.hh
+++ b/dune/stuff/functions/random_ellipsoids_function.hh
@@ -41,7 +41,7 @@ struct Ellipsoid
}
return DSC::FloatCmp::le(sum, 1.);
}
- bool intersects_cube(DomainType ll, DomainType ur) const
+ bool intersects_cube(DomainType /*ll*/, DomainType /*ur*/) const
{
DUNE_THROW(NotImplemented, "");
}
@@ -197,7 +197,7 @@ public:
EllipsoidType child = parent;
const double scale = std::pow(ellipsoid_cfg.get("ellipsoids.recursion_scale", 0.5), current_level);
const auto displace = [&](DomainFieldType& coord) {
- const auto disp = dist_rng() * DSC::sign(sign_rng());
+ const auto disp = dist_rng() * DSC::signum(sign_rng());
coord += disp;
DSC_LOG_DEBUG_0 << disp << ";";
};
diff --git a/dune/stuff/functions/spe10model2.hh b/dune/stuff/functions/spe10model2.hh
index e30a6b36a4dbdce36e64776bc3d79dd30ff41ed5..ad402bb07a08c14353f3d36250da4dbd4ad740d0 100644
--- a/dune/stuff/functions/spe10model2.hh
+++ b/dune/stuff/functions/spe10model2.hh
@@ -82,7 +82,7 @@ public:
}
}
- virtual size_t order() const
+ virtual size_t order() const override
{
return 0u;
}
diff --git a/dune/stuff/test/functions_ESV2007.cc b/dune/stuff/test/functions_ESV2007.cc
new file mode 100644
index 0000000000000000000000000000000000000000..2a447a38281c0d42a868a740e919344dec2b1e12
--- /dev/null
+++ b/dune/stuff/test/functions_ESV2007.cc
@@ -0,0 +1,78 @@
+// This file is part of the dune-stuff project:
+// https://github.com/wwu-numerik/dune-stuff
+// Copyright holders: Rene Milk, Felix Schindler
+// License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+
+#include "main.hxx"
+
+#include <memory>
+
+#include <dune/common/exceptions.hh>
+
+#include <dune/stuff/functions/interfaces.hh>
+#include <dune/stuff/functions/ESV2007.hh>
+
+
+// we need this nasty code generation because the testing::Types< ... > only accepts 50 arguments
+// and all combinations of functions and entities and dimensions and fieldtypes would be way too much
+#define TEST_STRUCT_GENERATOR(ftype, etype) \
+ template <class LocalizableFunctionType> \
+ struct ftype##etype##Test : public ::testing::Test \
+ { \
+ typedef typename LocalizableFunctionType::EntityType EntityType; \
+ typedef typename LocalizableFunctionType::LocalfunctionType LocalfunctionType; \
+ typedef typename LocalizableFunctionType::DomainFieldType DomainFieldType; \
+ static const size_t dimDomain = LocalizableFunctionType::dimDomain; \
+ typedef typename LocalizableFunctionType::DomainType DomainType; \
+ typedef typename LocalizableFunctionType::RangeFieldType RangeFieldType; \
+ static const size_t dimRange = LocalizableFunctionType::dimRange; \
+ static const size_t dimRangeCols = LocalizableFunctionType::dimRangeCols; \
+ typedef typename LocalizableFunctionType::RangeType RangeType; \
+ typedef typename LocalizableFunctionType::JacobianRangeType JacobianRangeType; \
+ \
+ void check() const \
+ { \
+ const std::unique_ptr<const LocalizableFunctionType> function( \
+ LocalizableFunctionType::create(LocalizableFunctionType::default_config())); \
+ } \
+ };
+// TEST_STRUCT_GENERATOR
+
+
+#if HAVE_DUNE_GRID
+
+#include <dune/grid/yaspgrid.hh>
+
+typedef Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<double, 2>>::Codim<0>::Entity DuneYaspGrid2dEntityType;
+
+typedef testing::
+ Types<Dune::Stuff::Functions::ESV2007::Testcase1Force<DuneYaspGrid2dEntityType, double, 2, double, 1, 1>,
+ Dune::Stuff::Functions::ESV2007::Testcase1ExactSolution<DuneYaspGrid2dEntityType, double, 2, double, 1, 1>>
+ ESV2007FunctionYaspGridEntityTypes;
+
+TEST_STRUCT_GENERATOR(CheckerboardFunction, YaspGridEntity)
+TYPED_TEST_CASE(CheckerboardFunctionYaspGridEntityTest, ESV2007FunctionYaspGridEntityTypes);
+TYPED_TEST(CheckerboardFunctionYaspGridEntityTest, provides_required_methods)
+{
+ this->check();
+}
+
+#if HAVE_ALUGRID
+#include <dune/grid/alugrid.hh>
+
+typedef Dune::ALUGrid<2, 2, Dune::simplex, Dune::nonconforming>::Codim<0>::Entity DuneAluSimplexGrid2dEntityType;
+
+typedef testing::
+ Types<Dune::Stuff::Functions::ESV2007::Testcase1Force<DuneAluSimplexGrid2dEntityType, double, 2, double, 1, 1>,
+ Dune::Stuff::Functions::ESV2007::Testcase1ExactSolution<DuneAluSimplexGrid2dEntityType, double, 2, double, 1,
+ 1>> ESV2007FunctionAluGridEntityTypes;
+
+TEST_STRUCT_GENERATOR(CheckerboardFunction, AluGridEntity)
+TYPED_TEST_CASE(CheckerboardFunctionAluGridEntityTest, ESV2007FunctionAluGridEntityTypes);
+TYPED_TEST(CheckerboardFunctionAluGridEntityTest, provides_required_methods)
+{
+ this->check();
+}
+
+#endif // HAVE_ALUGRID_SERIAL || HAVE_ALUGRID_PARALLEL
+#endif // HAVE_DUNE_GRID
diff --git a/dune/stuff/test/functions_checkerboard.cc b/dune/stuff/test/functions_checkerboard.cc
index e08ee24242c5df18acd31e9bb7606e5702398cb2..a8c923e9abdbb0d24f8a33e7651cc35d319bf2aa 100644
--- a/dune/stuff/test/functions_checkerboard.cc
+++ b/dune/stuff/test/functions_checkerboard.cc
@@ -41,54 +41,11 @@
#if HAVE_DUNE_GRID
-#include <dune/grid/sgrid.hh>
-
-typedef Dune::SGrid<1, 1>::Codim<0>::Entity DuneSGrid1dEntityType;
-typedef Dune::SGrid<2, 2>::Codim<0>::Entity DuneSGrid2dEntityType;
-typedef Dune::SGrid<3, 3>::Codim<0>::Entity DuneSGrid3dEntityType;
-
-typedef testing::Types<Dune::Stuff::Functions::Checkerboard<DuneSGrid1dEntityType, double, 1, double, 1, 1>,
- Dune::Stuff::Functions::Checkerboard<DuneSGrid1dEntityType, double, 1, double, 1, 2>,
- Dune::Stuff::Functions::Checkerboard<DuneSGrid1dEntityType, double, 1, double, 1, 3>,
- Dune::Stuff::Functions::Checkerboard<DuneSGrid1dEntityType, double, 1, double, 2, 1>,
- Dune::Stuff::Functions::Checkerboard<DuneSGrid1dEntityType, double, 1, double, 2, 2>,
- Dune::Stuff::Functions::Checkerboard<DuneSGrid1dEntityType, double, 1, double, 2, 3>,
- Dune::Stuff::Functions::Checkerboard<DuneSGrid1dEntityType, double, 1, double, 3, 1>,
- Dune::Stuff::Functions::Checkerboard<DuneSGrid1dEntityType, double, 1, double, 3, 2>,
- Dune::Stuff::Functions::Checkerboard<DuneSGrid1dEntityType, double, 1, double, 3, 3>,
- Dune::Stuff::Functions::Checkerboard<DuneSGrid2dEntityType, double, 2, double, 1, 1>,
- Dune::Stuff::Functions::Checkerboard<DuneSGrid2dEntityType, double, 2, double, 1, 2>,
- Dune::Stuff::Functions::Checkerboard<DuneSGrid2dEntityType, double, 2, double, 1, 3>,
- Dune::Stuff::Functions::Checkerboard<DuneSGrid2dEntityType, double, 2, double, 2, 1>,
- Dune::Stuff::Functions::Checkerboard<DuneSGrid2dEntityType, double, 2, double, 2, 2>,
- Dune::Stuff::Functions::Checkerboard<DuneSGrid2dEntityType, double, 2, double, 2, 3>,
- Dune::Stuff::Functions::Checkerboard<DuneSGrid2dEntityType, double, 2, double, 3, 1>,
- Dune::Stuff::Functions::Checkerboard<DuneSGrid2dEntityType, double, 2, double, 3, 2>,
- Dune::Stuff::Functions::Checkerboard<DuneSGrid2dEntityType, double, 2, double, 3, 3>,
- Dune::Stuff::Functions::Checkerboard<DuneSGrid3dEntityType, double, 3, double, 1, 1>,
- Dune::Stuff::Functions::Checkerboard<DuneSGrid3dEntityType, double, 3, double, 1, 2>,
- Dune::Stuff::Functions::Checkerboard<DuneSGrid3dEntityType, double, 3, double, 1, 3>,
- Dune::Stuff::Functions::Checkerboard<DuneSGrid3dEntityType, double, 3, double, 2, 1>,
- Dune::Stuff::Functions::Checkerboard<DuneSGrid3dEntityType, double, 3, double, 2, 2>,
- Dune::Stuff::Functions::Checkerboard<DuneSGrid3dEntityType, double, 3, double, 2, 3>,
- Dune::Stuff::Functions::Checkerboard<DuneSGrid3dEntityType, double, 3, double, 3, 1>,
- Dune::Stuff::Functions::Checkerboard<DuneSGrid3dEntityType, double, 3, double, 3, 2>,
- Dune::Stuff::Functions::Checkerboard<DuneSGrid3dEntityType, double, 3, double, 3, 3>,
- Dune::Stuff::Functions::Checkerboard<DuneSGrid3dEntityType, double, 3, std::complex<double>, 3,
- 3>> CheckerboardFunctionSGridEntityTypes;
-
-TEST_STRUCT_GENERATOR(CheckerboardFunction, SGridEntity)
-TYPED_TEST_CASE(CheckerboardFunctionSGridEntityTest, CheckerboardFunctionSGridEntityTypes);
-TYPED_TEST(CheckerboardFunctionSGridEntityTest, provides_required_methods)
-{
- this->check();
-}
-
#include <dune/grid/yaspgrid.hh>
-typedef Dune::YaspGrid<1>::Codim<0>::Entity DuneYaspGrid1dEntityType;
-typedef Dune::YaspGrid<2>::Codim<0>::Entity DuneYaspGrid2dEntityType;
-typedef Dune::YaspGrid<3>::Codim<0>::Entity DuneYaspGrid3dEntityType;
+typedef Dune::YaspGrid<1, Dune::EquidistantOffsetCoordinates<double, 1>>::Codim<0>::Entity DuneYaspGrid1dEntityType;
+typedef Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<double, 2>>::Codim<0>::Entity DuneYaspGrid2dEntityType;
+typedef Dune::YaspGrid<3, Dune::EquidistantOffsetCoordinates<double, 3>>::Codim<0>::Entity DuneYaspGrid3dEntityType;
typedef testing::Types<Dune::Stuff::Functions::Checkerboard<DuneYaspGrid1dEntityType, double, 1, double, 1, 1>,
Dune::Stuff::Functions::Checkerboard<DuneYaspGrid1dEntityType, double, 1, double, 1, 2>,
diff --git a/dune/stuff/test/functions_combined.cc b/dune/stuff/test/functions_combined.cc
index 2814075ff2c6f2363ce15ee909ffc757a0aca6a3..ed77358e65fdf69ad9c98d6bea7cb035dd3e16f4 100644
--- a/dune/stuff/test/functions_combined.cc
+++ b/dune/stuff/test/functions_combined.cc
@@ -10,7 +10,7 @@
#include <boost/numeric/conversion/cast.hpp>
#if HAVE_DUNE_GRID
-#include <dune/grid/sgrid.hh>
+#include <dune/grid/yaspgrid.hh>
#endif
#include <dune/geometry/quadraturerules.hh>
@@ -46,10 +46,12 @@ public:
template <class DimDomain>
class DifferenceFunctionTest
- : public FunctionTest<typename DifferenceFunctionType<SGrid<DimDomain::value, DimDomain::value>>::value>
+ : public FunctionTest<
+ typename DifferenceFunctionType<YaspGrid<DimDomain::value,
+ EquidistantOffsetCoordinates<double, DimDomain::value>>>::value>
{
protected:
- typedef SGrid<DimDomain::value, DimDomain::value> GridType;
+ typedef YaspGrid<DimDomain::value, EquidistantOffsetCoordinates<double, DimDomain::value>> GridType;
typedef typename DifferenceFunctionType<GridType>::value FunctionType;
static std::shared_ptr<GridType> create_grid()
@@ -59,8 +61,9 @@ protected:
static std::unique_ptr<FunctionType> create(const double ll, const double rr)
{
- typedef typename DifferenceFunctionType<SGrid<DimDomain::value, DimDomain::value>>::ConstantFunctionType
- ConstantFunctionType;
+ typedef typename DifferenceFunctionType<YaspGrid<DimDomain::value,
+ EquidistantOffsetCoordinates<double, DimDomain::value>>>::
+ ConstantFunctionType ConstantFunctionType;
auto left = std::make_shared<ConstantFunctionType>(ll);
auto right = std::make_shared<ConstantFunctionType>(rr);
return std::unique_ptr<FunctionType>(new FunctionType(left, right));
diff --git a/dune/stuff/test/functions_constant.cc b/dune/stuff/test/functions_constant.cc
index 8c5b21e97eabc0fda8d50a81c5d82e8144fba383..3271708153ea67371fdf69aa404823160d8d60d4 100644
--- a/dune/stuff/test/functions_constant.cc
+++ b/dune/stuff/test/functions_constant.cc
@@ -40,54 +40,11 @@
#if HAVE_DUNE_GRID
-#include <dune/grid/sgrid.hh>
-
-typedef Dune::SGrid<1, 1>::Codim<0>::Entity DuneSGrid1dEntityType;
-typedef Dune::SGrid<2, 2>::Codim<0>::Entity DuneSGrid2dEntityType;
-typedef Dune::SGrid<3, 3>::Codim<0>::Entity DuneSGrid3dEntityType;
-
-typedef testing::Types<Dune::Stuff::Functions::Constant<DuneSGrid1dEntityType, double, 1, double, 1, 1>,
- Dune::Stuff::Functions::Constant<DuneSGrid1dEntityType, double, 1, double, 1, 2>,
- Dune::Stuff::Functions::Constant<DuneSGrid1dEntityType, double, 1, double, 1, 3>,
- Dune::Stuff::Functions::Constant<DuneSGrid1dEntityType, double, 1, double, 2, 1>,
- Dune::Stuff::Functions::Constant<DuneSGrid1dEntityType, double, 1, double, 2, 2>,
- Dune::Stuff::Functions::Constant<DuneSGrid1dEntityType, double, 1, double, 2, 3>,
- Dune::Stuff::Functions::Constant<DuneSGrid1dEntityType, double, 1, double, 3, 1>,
- Dune::Stuff::Functions::Constant<DuneSGrid1dEntityType, double, 1, double, 3, 2>,
- Dune::Stuff::Functions::Constant<DuneSGrid1dEntityType, double, 1, double, 3, 3>,
- Dune::Stuff::Functions::Constant<DuneSGrid2dEntityType, double, 2, double, 1, 1>,
- Dune::Stuff::Functions::Constant<DuneSGrid2dEntityType, double, 2, double, 1, 2>,
- Dune::Stuff::Functions::Constant<DuneSGrid2dEntityType, double, 2, double, 1, 3>,
- Dune::Stuff::Functions::Constant<DuneSGrid2dEntityType, double, 2, double, 2, 1>,
- Dune::Stuff::Functions::Constant<DuneSGrid2dEntityType, double, 2, double, 2, 2>,
- Dune::Stuff::Functions::Constant<DuneSGrid2dEntityType, double, 2, double, 2, 3>,
- Dune::Stuff::Functions::Constant<DuneSGrid2dEntityType, double, 2, double, 3, 1>,
- Dune::Stuff::Functions::Constant<DuneSGrid2dEntityType, double, 2, double, 3, 2>,
- Dune::Stuff::Functions::Constant<DuneSGrid2dEntityType, double, 2, double, 3, 3>,
- Dune::Stuff::Functions::Constant<DuneSGrid3dEntityType, double, 3, double, 1, 1>,
- Dune::Stuff::Functions::Constant<DuneSGrid3dEntityType, double, 3, double, 1, 2>,
- Dune::Stuff::Functions::Constant<DuneSGrid3dEntityType, double, 3, double, 1, 3>,
- Dune::Stuff::Functions::Constant<DuneSGrid3dEntityType, double, 3, double, 2, 1>,
- Dune::Stuff::Functions::Constant<DuneSGrid3dEntityType, double, 3, double, 2, 2>,
- Dune::Stuff::Functions::Constant<DuneSGrid3dEntityType, double, 3, double, 2, 3>,
- Dune::Stuff::Functions::Constant<DuneSGrid3dEntityType, double, 3, double, 3, 1>,
- Dune::Stuff::Functions::Constant<DuneSGrid3dEntityType, double, 3, double, 3, 2>,
- Dune::Stuff::Functions::Constant<DuneSGrid3dEntityType, double, 3, double, 3, 3>,
- Dune::Stuff::Functions::Constant<DuneSGrid3dEntityType, double, 3, std::complex<double>, 3, 3>>
- ConstantFunctionSGridEntityTypes;
-
-TEST_STRUCT_GENERATOR(ConstantFunction, SGridEntity)
-TYPED_TEST_CASE(ConstantFunctionSGridEntityTest, ConstantFunctionSGridEntityTypes);
-TYPED_TEST(ConstantFunctionSGridEntityTest, provides_required_methods)
-{
- this->check();
-}
-
#include <dune/grid/yaspgrid.hh>
-typedef Dune::YaspGrid<1>::Codim<0>::Entity DuneYaspGrid1dEntityType;
-typedef Dune::YaspGrid<2>::Codim<0>::Entity DuneYaspGrid2dEntityType;
-typedef Dune::YaspGrid<3>::Codim<0>::Entity DuneYaspGrid3dEntityType;
+typedef Dune::YaspGrid<1, Dune::EquidistantOffsetCoordinates<double, 1>>::Codim<0>::Entity DuneYaspGrid1dEntityType;
+typedef Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<double, 2>>::Codim<0>::Entity DuneYaspGrid2dEntityType;
+typedef Dune::YaspGrid<3, Dune::EquidistantOffsetCoordinates<double, 3>>::Codim<0>::Entity DuneYaspGrid3dEntityType;
typedef testing::Types<Dune::Stuff::Functions::Constant<DuneYaspGrid1dEntityType, double, 1, double, 1, 1>,
Dune::Stuff::Functions::Constant<DuneYaspGrid1dEntityType, double, 1, double, 1, 2>,
diff --git a/dune/stuff/test/functions_expression.cc b/dune/stuff/test/functions_expression.cc
index ab8e93aa4fdc0577e3948d5a04036cfb404bd9fe..cb9dfe26fa3299d1b1ca7b86f7a77bcb9fa6b73d 100644
--- a/dune/stuff/test/functions_expression.cc
+++ b/dune/stuff/test/functions_expression.cc
@@ -48,54 +48,11 @@
#if HAVE_DUNE_GRID
-#include <dune/grid/sgrid.hh>
-
-typedef Dune::SGrid<1, 1>::Codim<0>::Entity DuneSGrid1dEntityType;
-typedef Dune::SGrid<2, 2>::Codim<0>::Entity DuneSGrid2dEntityType;
-typedef Dune::SGrid<3, 3>::Codim<0>::Entity DuneSGrid3dEntityType;
-
-// the matrix valued version is missing the create() and default_config() method
-typedef testing::Types<Dune::Stuff::Functions::Expression<DuneSGrid1dEntityType, double, 1, double, 1, 1>,
- Dune::Stuff::Functions::Expression<DuneSGrid1dEntityType, double, 1, double, 1, 2>,
- Dune::Stuff::Functions::Expression<DuneSGrid1dEntityType, double, 1, double, 1, 3>,
- Dune::Stuff::Functions::Expression<DuneSGrid1dEntityType, double, 1, double, 2, 1>,
- Dune::Stuff::Functions::Expression<DuneSGrid1dEntityType, double, 1, double, 2, 2>,
- Dune::Stuff::Functions::Expression<DuneSGrid1dEntityType, double, 1, double, 2, 3>,
- Dune::Stuff::Functions::Expression<DuneSGrid1dEntityType, double, 1, double, 3, 1>,
- Dune::Stuff::Functions::Expression<DuneSGrid1dEntityType, double, 1, double, 3, 2>,
- Dune::Stuff::Functions::Expression<DuneSGrid1dEntityType, double, 1, double, 3, 3>,
- Dune::Stuff::Functions::Expression<DuneSGrid2dEntityType, double, 2, double, 1, 1>,
- Dune::Stuff::Functions::Expression<DuneSGrid2dEntityType, double, 2, double, 1, 2>,
- Dune::Stuff::Functions::Expression<DuneSGrid2dEntityType, double, 2, double, 1, 3>,
- Dune::Stuff::Functions::Expression<DuneSGrid2dEntityType, double, 2, double, 2, 1>,
- Dune::Stuff::Functions::Expression<DuneSGrid2dEntityType, double, 2, double, 2, 2>,
- Dune::Stuff::Functions::Expression<DuneSGrid2dEntityType, double, 2, double, 2, 3>,
- Dune::Stuff::Functions::Expression<DuneSGrid2dEntityType, double, 2, double, 3, 1>,
- Dune::Stuff::Functions::Expression<DuneSGrid2dEntityType, double, 2, double, 3, 2>,
- Dune::Stuff::Functions::Expression<DuneSGrid2dEntityType, double, 2, double, 3, 3>,
- Dune::Stuff::Functions::Expression<DuneSGrid3dEntityType, double, 3, double, 1, 1>,
- Dune::Stuff::Functions::Expression<DuneSGrid3dEntityType, double, 3, double, 1, 2>,
- Dune::Stuff::Functions::Expression<DuneSGrid3dEntityType, double, 3, double, 1, 3>,
- Dune::Stuff::Functions::Expression<DuneSGrid3dEntityType, double, 3, double, 2, 1>,
- Dune::Stuff::Functions::Expression<DuneSGrid3dEntityType, double, 3, double, 2, 2>,
- Dune::Stuff::Functions::Expression<DuneSGrid3dEntityType, double, 3, double, 2, 3>,
- Dune::Stuff::Functions::Expression<DuneSGrid3dEntityType, double, 3, double, 3, 1>,
- Dune::Stuff::Functions::Expression<DuneSGrid3dEntityType, double, 3, double, 3, 2>,
- Dune::Stuff::Functions::Expression<DuneSGrid3dEntityType, double, 3, double, 3, 3>>
- ExpressionFunctionSGridEntityTypes;
-
-TEST_STRUCT_GENERATOR(ExpressionFunction, SGridEntity)
-TYPED_TEST_CASE(ExpressionFunctionSGridEntityTest, ExpressionFunctionSGridEntityTypes);
-TYPED_TEST(ExpressionFunctionSGridEntityTest, provides_required_methods)
-{
- this->check();
-}
-
#include <dune/grid/yaspgrid.hh>
-typedef Dune::YaspGrid<1>::Codim<0>::Entity DuneYaspGrid1dEntityType;
-typedef Dune::YaspGrid<2>::Codim<0>::Entity DuneYaspGrid2dEntityType;
-typedef Dune::YaspGrid<3>::Codim<0>::Entity DuneYaspGrid3dEntityType;
+typedef Dune::YaspGrid<1, Dune::EquidistantOffsetCoordinates<double, 1>>::Codim<0>::Entity DuneYaspGrid1dEntityType;
+typedef Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<double, 2>>::Codim<0>::Entity DuneYaspGrid2dEntityType;
+typedef Dune::YaspGrid<3, Dune::EquidistantOffsetCoordinates<double, 3>>::Codim<0>::Entity DuneYaspGrid3dEntityType;
typedef testing::Types<Dune::Stuff::Functions::Expression<DuneYaspGrid1dEntityType, double, 1, double, 1, 1>,
Dune::Stuff::Functions::Expression<DuneYaspGrid1dEntityType, double, 1, double, 1, 2>,
diff --git a/dune/stuff/test/functions_flattop.cc b/dune/stuff/test/functions_flattop.cc
index 36d9529a6546c571a638c6ea75f88683223535b8..ad5389f507973424ace4659b5f0773c6161cd976 100644
--- a/dune/stuff/test/functions_flattop.cc
+++ b/dune/stuff/test/functions_flattop.cc
@@ -14,7 +14,7 @@
#include <boost/numeric/conversion/cast.hpp>
#if HAVE_DUNE_GRID
-#include <dune/grid/sgrid.hh>
+#include <dune/grid/yaspgrid.hh>
#endif
#include <dune/geometry/quadraturerules.hh>
@@ -50,10 +50,12 @@ public:
template <class DimDomain>
class FlatTopFunctionTest
- : public FunctionTest<typename FlatTopFunctionType<SGrid<DimDomain::value, DimDomain::value>>::value>
+ : public FunctionTest<
+ typename FlatTopFunctionType<YaspGrid<DimDomain::value,
+ EquidistantOffsetCoordinates<double, DimDomain::value>>>::value>
{
protected:
- typedef SGrid<DimDomain::value, DimDomain::value> GridType;
+ typedef YaspGrid<DimDomain::value, EquidistantOffsetCoordinates<double, DimDomain::value>> GridType;
typedef typename FlatTopFunctionType<GridType>::value FunctionType;
static std::shared_ptr<GridType> create_grid()
diff --git a/dune/stuff/test/functions_functions.cc b/dune/stuff/test/functions_functions.cc
index ab2a565973831ef0bf8df25c92532a2619a6dd2d..9b675548656045abb9c97cc891662e344e43c44c 100644
--- a/dune/stuff/test/functions_functions.cc
+++ b/dune/stuff/test/functions_functions.cc
@@ -45,54 +45,11 @@
#if HAVE_DUNE_GRID
-#include <dune/grid/sgrid.hh>
-
-typedef Dune::SGrid<1, 1>::Codim<0>::Entity DuneSGrid1dEntityType;
-typedef Dune::SGrid<2, 2>::Codim<0>::Entity DuneSGrid2dEntityType;
-typedef Dune::SGrid<3, 3>::Codim<0>::Entity DuneSGrid3dEntityType;
-
-typedef testing::Types<Dune::Stuff::Functions::Constant<DuneSGrid1dEntityType, double, 1, double, 1, 1>,
- Dune::Stuff::Functions::Constant<DuneSGrid1dEntityType, double, 1, double, 1, 2>,
- Dune::Stuff::Functions::Constant<DuneSGrid1dEntityType, double, 1, double, 1, 3>,
- Dune::Stuff::Functions::Constant<DuneSGrid1dEntityType, double, 1, double, 2, 1>,
- Dune::Stuff::Functions::Constant<DuneSGrid1dEntityType, double, 1, double, 2, 2>,
- Dune::Stuff::Functions::Constant<DuneSGrid1dEntityType, double, 1, double, 2, 3>,
- Dune::Stuff::Functions::Constant<DuneSGrid1dEntityType, double, 1, double, 3, 1>,
- Dune::Stuff::Functions::Constant<DuneSGrid1dEntityType, double, 1, double, 3, 2>,
- Dune::Stuff::Functions::Constant<DuneSGrid1dEntityType, double, 1, double, 3, 3>,
- Dune::Stuff::Functions::Constant<DuneSGrid2dEntityType, double, 2, double, 1, 1>,
- Dune::Stuff::Functions::Constant<DuneSGrid2dEntityType, double, 2, double, 1, 2>,
- Dune::Stuff::Functions::Constant<DuneSGrid2dEntityType, double, 2, double, 1, 3>,
- Dune::Stuff::Functions::Constant<DuneSGrid2dEntityType, double, 2, double, 2, 1>,
- Dune::Stuff::Functions::Constant<DuneSGrid2dEntityType, double, 2, double, 2, 2>,
- Dune::Stuff::Functions::Constant<DuneSGrid2dEntityType, double, 2, double, 2, 3>,
- Dune::Stuff::Functions::Constant<DuneSGrid2dEntityType, double, 2, double, 3, 1>,
- Dune::Stuff::Functions::Constant<DuneSGrid2dEntityType, double, 2, double, 3, 2>,
- Dune::Stuff::Functions::Constant<DuneSGrid2dEntityType, double, 2, double, 3, 3>,
- Dune::Stuff::Functions::Constant<DuneSGrid3dEntityType, double, 3, double, 1, 1>,
- Dune::Stuff::Functions::Constant<DuneSGrid3dEntityType, double, 3, double, 1, 2>,
- Dune::Stuff::Functions::Constant<DuneSGrid3dEntityType, double, 3, double, 1, 3>,
- Dune::Stuff::Functions::Constant<DuneSGrid3dEntityType, double, 3, double, 2, 1>,
- Dune::Stuff::Functions::Constant<DuneSGrid3dEntityType, double, 3, double, 2, 2>,
- Dune::Stuff::Functions::Constant<DuneSGrid3dEntityType, double, 3, double, 2, 3>,
- Dune::Stuff::Functions::Constant<DuneSGrid3dEntityType, double, 3, double, 3, 1>,
- Dune::Stuff::Functions::Constant<DuneSGrid3dEntityType, double, 3, double, 3, 2>,
- Dune::Stuff::Functions::Constant<DuneSGrid3dEntityType, double, 3, double, 3, 3>,
- Dune::Stuff::Functions::Constant<DuneSGrid3dEntityType, double, 3, std::complex<double>, 3, 3>>
- FunctionsSGridEntityTypes;
-
-TEST_STRUCT_GENERATOR(Functions, SGridEntity)
-TYPED_TEST_CASE(FunctionsSGridEntityTest, FunctionsSGridEntityTypes);
-TYPED_TEST(FunctionsSGridEntityTest, provides_required_methods)
-{
- this->check();
-}
-
#include <dune/grid/yaspgrid.hh>
-typedef Dune::YaspGrid<1>::Codim<0>::Entity DuneYaspGrid1dEntityType;
-typedef Dune::YaspGrid<2>::Codim<0>::Entity DuneYaspGrid2dEntityType;
-typedef Dune::YaspGrid<3>::Codim<0>::Entity DuneYaspGrid3dEntityType;
+typedef Dune::YaspGrid<1, Dune::EquidistantOffsetCoordinates<double, 1>>::Codim<0>::Entity DuneYaspGrid1dEntityType;
+typedef Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<double, 2>>::Codim<0>::Entity DuneYaspGrid2dEntityType;
+typedef Dune::YaspGrid<3, Dune::EquidistantOffsetCoordinates<double, 3>>::Codim<0>::Entity DuneYaspGrid3dEntityType;
typedef testing::Types<Dune::Stuff::Functions::Constant<DuneYaspGrid1dEntityType, double, 1, double, 1, 1>,
Dune::Stuff::Functions::Constant<DuneYaspGrid1dEntityType, double, 1, double, 1, 2>,
diff --git a/dune/stuff/test/functions_global.cc b/dune/stuff/test/functions_global.cc
index c32997cab45590fa0c960d546867c39e39345445..b4ccfa5c026fe95065be8989e4431081be577e4a 100644
--- a/dune/stuff/test/functions_global.cc
+++ b/dune/stuff/test/functions_global.cc
@@ -45,54 +45,11 @@
#if HAVE_DUNE_GRID
-#include <dune/grid/sgrid.hh>
-
-typedef Dune::SGrid<1, 1>::Codim<0>::Entity DuneSGrid1dEntityType;
-typedef Dune::SGrid<2, 2>::Codim<0>::Entity DuneSGrid2dEntityType;
-typedef Dune::SGrid<3, 3>::Codim<0>::Entity DuneSGrid3dEntityType;
-
-typedef testing::Types<Dune::Stuff::GlobalLambdaFunction<DuneSGrid1dEntityType, double, 1, double, 1, 1>,
- Dune::Stuff::GlobalLambdaFunction<DuneSGrid1dEntityType, double, 1, double, 1, 2>,
- Dune::Stuff::GlobalLambdaFunction<DuneSGrid1dEntityType, double, 1, double, 1, 3>,
- Dune::Stuff::GlobalLambdaFunction<DuneSGrid1dEntityType, double, 1, double, 2, 1>,
- Dune::Stuff::GlobalLambdaFunction<DuneSGrid1dEntityType, double, 1, double, 2, 2>,
- Dune::Stuff::GlobalLambdaFunction<DuneSGrid1dEntityType, double, 1, double, 2, 3>,
- Dune::Stuff::GlobalLambdaFunction<DuneSGrid1dEntityType, double, 1, double, 3, 1>,
- Dune::Stuff::GlobalLambdaFunction<DuneSGrid1dEntityType, double, 1, double, 3, 2>,
- Dune::Stuff::GlobalLambdaFunction<DuneSGrid1dEntityType, double, 1, double, 3, 3>,
- Dune::Stuff::GlobalLambdaFunction<DuneSGrid2dEntityType, double, 2, double, 1, 1>,
- Dune::Stuff::GlobalLambdaFunction<DuneSGrid2dEntityType, double, 2, double, 1, 2>,
- Dune::Stuff::GlobalLambdaFunction<DuneSGrid2dEntityType, double, 2, double, 1, 3>,
- Dune::Stuff::GlobalLambdaFunction<DuneSGrid2dEntityType, double, 2, double, 2, 1>,
- Dune::Stuff::GlobalLambdaFunction<DuneSGrid2dEntityType, double, 2, double, 2, 2>,
- Dune::Stuff::GlobalLambdaFunction<DuneSGrid2dEntityType, double, 2, double, 2, 3>,
- Dune::Stuff::GlobalLambdaFunction<DuneSGrid2dEntityType, double, 2, double, 3, 1>,
- Dune::Stuff::GlobalLambdaFunction<DuneSGrid2dEntityType, double, 2, double, 3, 2>,
- Dune::Stuff::GlobalLambdaFunction<DuneSGrid2dEntityType, double, 2, double, 3, 3>,
- Dune::Stuff::GlobalLambdaFunction<DuneSGrid3dEntityType, double, 3, double, 1, 1>,
- Dune::Stuff::GlobalLambdaFunction<DuneSGrid3dEntityType, double, 3, double, 1, 2>,
- Dune::Stuff::GlobalLambdaFunction<DuneSGrid3dEntityType, double, 3, double, 1, 3>,
- Dune::Stuff::GlobalLambdaFunction<DuneSGrid3dEntityType, double, 3, double, 2, 1>,
- Dune::Stuff::GlobalLambdaFunction<DuneSGrid3dEntityType, double, 3, double, 2, 2>,
- Dune::Stuff::GlobalLambdaFunction<DuneSGrid3dEntityType, double, 3, double, 2, 3>,
- Dune::Stuff::GlobalLambdaFunction<DuneSGrid3dEntityType, double, 3, double, 3, 1>,
- Dune::Stuff::GlobalLambdaFunction<DuneSGrid3dEntityType, double, 3, double, 3, 2>,
- Dune::Stuff::GlobalLambdaFunction<DuneSGrid3dEntityType, double, 3, double, 3, 3>,
- Dune::Stuff::GlobalLambdaFunction<DuneSGrid3dEntityType, double, 3, std::complex<double>, 3, 3>>
- GlobalLambdaFunctionSGridEntityTypes;
-
-TEST_STRUCT_GENERATOR(GlobalLambdaFunctionTest, SGridEntity)
-TYPED_TEST_CASE(GlobalLambdaFunctionTestSGridEntityTest, GlobalLambdaFunctionSGridEntityTypes);
-TYPED_TEST(GlobalLambdaFunctionTestSGridEntityTest, provides_required_methods)
-{
- this->check();
-}
-
#include <dune/grid/yaspgrid.hh>
-typedef Dune::YaspGrid<1>::Codim<0>::Entity DuneYaspGrid1dEntityType;
-typedef Dune::YaspGrid<2>::Codim<0>::Entity DuneYaspGrid2dEntityType;
-typedef Dune::YaspGrid<3>::Codim<0>::Entity DuneYaspGrid3dEntityType;
+typedef Dune::YaspGrid<1, Dune::EquidistantOffsetCoordinates<double, 1>>::Codim<0>::Entity DuneYaspGrid1dEntityType;
+typedef Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<double, 2>>::Codim<0>::Entity DuneYaspGrid2dEntityType;
+typedef Dune::YaspGrid<3, Dune::EquidistantOffsetCoordinates<double, 3>>::Codim<0>::Entity DuneYaspGrid3dEntityType;
typedef testing::Types<Dune::Stuff::GlobalLambdaFunction<DuneYaspGrid1dEntityType, double, 1, double, 1, 1>,
Dune::Stuff::GlobalLambdaFunction<DuneYaspGrid1dEntityType, double, 1, double, 1, 2>,
diff --git a/dune/stuff/test/functions_indicator.cc b/dune/stuff/test/functions_indicator.cc
new file mode 100644
index 0000000000000000000000000000000000000000..c3a051abfe3b9834325582ed1a1b4f19931ba4a4
--- /dev/null
+++ b/dune/stuff/test/functions_indicator.cc
@@ -0,0 +1,82 @@
+// This file is part of the dune-stuff project:
+// https://github.com/wwu-numerik/dune-stuff
+// Copyright holders: Rene Milk, Felix Schindler
+// License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+
+#include "main.hxx"
+
+#include <memory>
+
+#include <dune/common/exceptions.hh>
+
+#include <dune/stuff/functions/interfaces.hh>
+#include <dune/stuff/functions/indicator.hh>
+
+
+// we need this nasty code generation because the testing::Types< ... > only accepts 50 arguments
+// and all combinations of functions and entities and dimensions and fieldtypes would be way too much
+#define TEST_STRUCT_GENERATOR(ftype, etype) \
+ template <class LocalizableFunctionType> \
+ struct ftype##etype##Test : public ::testing::Test \
+ { \
+ typedef typename LocalizableFunctionType::EntityType EntityType; \
+ typedef typename LocalizableFunctionType::LocalfunctionType LocalfunctionType; \
+ typedef typename LocalizableFunctionType::DomainFieldType DomainFieldType; \
+ static const size_t dimDomain = LocalizableFunctionType::dimDomain; \
+ typedef typename LocalizableFunctionType::DomainType DomainType; \
+ typedef typename LocalizableFunctionType::RangeFieldType RangeFieldType; \
+ static const size_t dimRange = LocalizableFunctionType::dimRange; \
+ static const size_t dimRangeCols = LocalizableFunctionType::dimRangeCols; \
+ typedef typename LocalizableFunctionType::RangeType RangeType; \
+ typedef typename LocalizableFunctionType::JacobianRangeType JacobianRangeType; \
+ \
+ void check() const \
+ { \
+ const std::unique_ptr<const LocalizableFunctionType> function( \
+ LocalizableFunctionType::create(LocalizableFunctionType::default_config())); \
+ } \
+ };
+// TEST_STRUCT_GENERATOR
+
+
+#if HAVE_DUNE_GRID
+
+#include <dune/grid/yaspgrid.hh>
+
+typedef Dune::YaspGrid<1, Dune::EquidistantOffsetCoordinates<double, 1>>::Codim<0>::Entity DuneYaspGrid1dEntityType;
+typedef Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<double, 2>>::Codim<0>::Entity DuneYaspGrid2dEntityType;
+typedef Dune::YaspGrid<3, Dune::EquidistantOffsetCoordinates<double, 3>>::Codim<0>::Entity DuneYaspGrid3dEntityType;
+
+typedef testing::Types<Dune::Stuff::Functions::Indicator<DuneYaspGrid1dEntityType, double, 1, double, 1, 1>,
+ Dune::Stuff::Functions::Indicator<DuneYaspGrid2dEntityType, double, 2, double, 1, 1>,
+ Dune::Stuff::Functions::Indicator<DuneYaspGrid3dEntityType, double, 3, double, 1, 1>>
+ IndicatorFunctionYaspGridEntityTypes;
+
+TEST_STRUCT_GENERATOR(IndicatorFunction, YaspGridEntity)
+TYPED_TEST_CASE(IndicatorFunctionYaspGridEntityTest, IndicatorFunctionYaspGridEntityTypes);
+TYPED_TEST(IndicatorFunctionYaspGridEntityTest, provides_required_methods)
+{
+ this->check();
+}
+
+#if HAVE_ALUGRID
+#include <dune/grid/alugrid.hh>
+
+typedef Dune::ALUGrid<2, 2, Dune::simplex, Dune::nonconforming>::Codim<0>::Entity DuneAluSimplexGrid2dEntityType;
+typedef Dune::ALUGrid<3, 3, Dune::simplex, Dune::nonconforming>::Codim<0>::Entity DuneAluSimplexGrid3dEntityType;
+typedef Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming>::Codim<0>::Entity DuneAluCubeGrid3dEntityType;
+
+typedef testing::Types<Dune::Stuff::Functions::Indicator<DuneAluSimplexGrid2dEntityType, double, 2, double, 1, 1>,
+ Dune::Stuff::Functions::Indicator<DuneAluSimplexGrid3dEntityType, double, 3, double, 1, 1>,
+ Dune::Stuff::Functions::Indicator<DuneAluCubeGrid3dEntityType, double, 3, double, 1, 1>>
+ IndicatorFunctionAluGridEntityTypes;
+
+TEST_STRUCT_GENERATOR(IndicatorFunction, AluGridEntity)
+TYPED_TEST_CASE(IndicatorFunctionAluGridEntityTest, IndicatorFunctionAluGridEntityTypes);
+TYPED_TEST(IndicatorFunctionAluGridEntityTest, provides_required_methods)
+{
+ this->check();
+}
+
+#endif // HAVE_ALUGRID_SERIAL || HAVE_ALUGRID_PARALLEL
+#endif // HAVE_DUNE_GRID
diff --git a/dune/stuff/test/functions_spe10.cc b/dune/stuff/test/functions_spe10.cc
index d70f56740078d326c0b7f15a6f0c2e980a479a5c..818a6096c3d5a87b12b5f94a8fa387480db41949 100644
--- a/dune/stuff/test/functions_spe10.cc
+++ b/dune/stuff/test/functions_spe10.cc
@@ -40,7 +40,7 @@
#if HAVE_DUNE_GRID
-//# include <dune/grid/sgrid.hh>
+//# include <dune/grid/yaspgrid.hh>
// typedef Dune::SGrid< 1, 1 >::Codim< 0 >::Entity DuneSGrid1dEntityType;
// typedef Dune::SGrid< 2, 2 >::Codim< 0 >::Entity DuneSGrid2dEntityType;
@@ -83,7 +83,7 @@
#include <dune/grid/yaspgrid.hh>
-typedef Dune::YaspGrid<3>::Codim<0>::Entity DuneYaspGrid3dEntityType;
+typedef Dune::YaspGrid<3, Dune::EquidistantOffsetCoordinates<double, 3>>::Codim<0>::Entity DuneYaspGrid3dEntityType;
typedef testing::Types<Dune::Stuff::Functions::Spe10::Model2<DuneYaspGrid3dEntityType, double, 3, double, 3, 3>>
Spe10Model2FunctionYaspGridEntityTypes;