FreeFEM3D (aka ff3d): MultiLinearExpression.cpp Source File

00001 //  This file is part of ff3d - http://www.freefem.org/ff3d
00002 //  Copyright (C) 2007 St�phane Del Pino
00003 
00004 //  This program is free software; you can redistribute it and/or modify
00005 //  it under the terms of the GNU General Public License as published by
00006 //  the Free Software Foundation; either version 2, or (at your option)
00007 //  any later version.
00008 
00009 //  This program is distributed in the hope that it will be useful,
00010 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00011 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012 //  GNU General Public License for more details.
00013 
00014 //  You should have received a copy of the GNU General Public License
00015 //  along with this program; if not, write to the Free Software Foundation,
00016 //  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  
00017 
00018 //  $Id: MultiLinearExpression.cpp,v 1.1 2007/06/10 14:59:57 delpinux Exp $
00019 
00020 #include <MultiLinearExpression.hpp>
00021 
00022 #include <EmbededFunctions.hpp>
00023 #include <Stringify.hpp>
00024 #include <ErrorHandler.hpp>
00025 
00026 #include <FunctionExpressionConstant.hpp>
00027 #include <FunctionExpressionBinaryOperation.hpp>
00028 
00029 std::ostream&
00030 MultiLinearExpression::
00031 put(std::ostream& os) const
00032 {
00033   bool notFirst = false;
00034   for (LinearListType::const_iterator i = __linearList.begin();
00035        i != __linearList.end(); ++i) {
00036     if (notFirst) {
00037       os << '*';
00038     } else {
00039       notFirst = true;
00040     }
00041     os << (*(*i));
00042   }
00043 
00044   for (FunctionListType::const_iterator i = __functionList.begin();
00045        i != __functionList.end(); ++i) {
00046     if (notFirst) {
00047       os << '*';
00048     } else {
00049       notFirst = true;
00050     }
00051     os << (*(*i));
00052   }
00053 
00054   for (RealListType::const_iterator i = __realList.begin();
00055        i != __realList.end(); ++i) {
00056     if (notFirst) {
00057       os << '*';
00058     } else {
00059       notFirst = true;
00060     }
00061     os << (*(*i));
00062   }
00063 
00064   return os;
00065 }
00066 
00067 
00068 ReferenceCounting<FunctionExpression>
00069 MultiLinearExpression::
00070 getFunction()
00071 {
00072   ReferenceCounting<FunctionExpression> f = 0;
00073   ReferenceCounting<RealExpression> r = 0;
00074 
00075   if (__realList.size() > 0) {
00076     RealListType::iterator i = __realList.begin();
00077     r = *i;
00078     i++;
00079     for ( ; i != __realList.end(); ++i) {
00080       r = new RealExpressionBinaryOperator<product>(r,*i);
00081     }
00082   }
00083   if (__functionList.size() > 0) {
00084     FunctionListType::iterator i = __functionList.begin();
00085     f = *i;
00086     ++i;
00087     for ( ; i != __functionList.end(); ++i) {
00088       f = new FunctionExpressionBinaryOperation(BinaryOperation::product, f, *i);
00089     }
00090   }
00091 
00092   for (LinearListType::iterator i = __linearList.begin();
00093        i != __linearList.end(); ++i) {
00094     switch ((*(*i)).formType()) {
00095     case (LinearExpression::elementary): {
00096       break;                    // nothing to do
00097     }
00098     case (LinearExpression::elementaryTimesFunction): {
00099       LinearExpressionElementaryTimesFunction& L
00100         = dynamic_cast<LinearExpressionElementaryTimesFunction&>(*(*i));
00101       if (f == 0) {
00102         f = L.function();
00103       } else {
00104         f = new FunctionExpressionBinaryOperation(BinaryOperation::product, f, L.function());
00105       }
00106       break;
00107     }
00108     case (LinearExpression::elementaryTimesReal): {
00109       LinearExpressionElementaryTimesReal& L
00110         = dynamic_cast<LinearExpressionElementaryTimesReal&>(*(*i));
00111       if (r == 0) {
00112         r = L.real_t();
00113       } else {
00114         r = new RealExpressionBinaryOperator<product>(r,L.real_t());
00115       }
00116       break;
00117     }
00118     case (LinearExpression::elementaryTimesFunctionOperator): {
00119       throw ErrorHandler(__FILE__,__LINE__,
00120                          "not implemented",
00121                          ErrorHandler::unexpected);
00122       break;
00123     }
00124     default: {
00125       throw ErrorHandler(__FILE__,__LINE__,
00126                          "not implemented",
00127                          ErrorHandler::unexpected);
00128     }
00129     }
00130   }
00131 
00132   if ((r == 0)&&(f == 0)) {
00133     r = new RealExpressionValue(1);
00134     f = new FunctionExpressionConstant(r);
00135   } else if (r == 0) {
00136     ; // everything is ok !
00137   } else if (f == 0){
00138     f = new FunctionExpressionConstant(r);
00139   } else {
00140     f = new FunctionExpressionBinaryOperation(BinaryOperation::product, f, new FunctionExpressionConstant(r));
00141   }
00142 
00143   return f;
00144 }
00145 
00146 const std::string&
00147 MultiLinearExpression::
00148 getUnknownName() const
00149 {
00150   LinearListType::const_iterator i = __linearList.begin();
00151   while (i != __linearList.end()) {
00152     const IntegratedExpression& IE =
00153       *(*i)->integrated()->integratedExpression();
00154 
00155     if (IE.integratedExpressionType()
00156         == IntegratedExpression::unknownFunction) {
00157       const IntegratedExpressionUnknown& I
00158         = dynamic_cast<const IntegratedExpressionUnknown&>(IE);
00159 
00160       return I.unknownName();
00161     }
00162     ++i;
00163   }
00164     
00165   throw ErrorHandler(__FILE__,__LINE__,
00166                      "Could not find unknown name",
00167                      ErrorHandler::unexpected);
00168   static std::string badName("unknown function");
00169   return badName;
00170 }
00171 
00172 
00173 const std::string&
00174 MultiLinearExpression::
00175 getTestFunctionName() const
00176 {
00177   ConstReferenceCounting<TestFunctionVariable> t = 0;
00178   LinearListType::const_iterator i = __linearList.begin();
00179   while (i != __linearList.end()) {
00180     const IntegratedExpression& IE =
00181       *(*i)->integrated()->integratedExpression();
00182 
00183     if (IE.integratedExpressionType()
00184         == IntegratedExpression::testFunction) {
00185       const IntegratedExpressionTest& I
00186         = dynamic_cast<const IntegratedExpressionTest&>(IE);
00187       return I.testFunctionName();
00188     }
00189     ++i;
00190   }
00191   throw ErrorHandler(__FILE__,__LINE__,
00192                      "Could not find testfunction name",
00193                      ErrorHandler::unexpected);
00194   static std::string badName("unknown function");
00195   return badName;
00196 }
00197 
00198 
00199 IntegratedOperatorExpression::OperatorType
00200 MultiLinearExpression::
00201 getUnknownOperator() const
00202 {
00203   IntegratedOperatorExpression::OperatorType t
00204     = IntegratedOperatorExpression::undefined;
00205 
00206   LinearListType::const_iterator i = __linearList.begin();
00207   while (i != __linearList.end()) {
00208     const IntegratedOperatorExpression& I = *(*i)->integrated();
00209 
00210     const IntegratedExpression& IE = *I.integratedExpression();
00211 
00212     if (IE.integratedExpressionType()
00213         == IntegratedExpression::unknownFunction) {
00214       t = I.operatorType();
00215     }
00216     ++i;
00217   }
00218 
00219   return t;
00220 }
00221 
00222 IntegratedOperatorExpression::OperatorType
00223 MultiLinearExpression::
00224 getTestFunctionOperator() const
00225 {
00226   IntegratedOperatorExpression::OperatorType t
00227     = IntegratedOperatorExpression::undefined;
00228 
00229   LinearListType::const_iterator i = __linearList.begin();
00230   while (i != __linearList.end()) {
00231     const IntegratedOperatorExpression& I = *(*i)->integrated();
00232 
00233     const IntegratedExpression& IE = *I.integratedExpression();
00234 
00235     if (IE.integratedExpressionType()
00236         == IntegratedExpression::testFunction) {
00237       t = I.operatorType();
00238     }
00239     ++i;
00240   }
00241 
00242   return t;
00243 }
00244 
00245 IntegratedOperatorExpression::OperatorType
00246 MultiLinearExpression::
00247 getFunctionOperator() const
00248 {
00249   IntegratedOperatorExpression::OperatorType t
00250     = IntegratedOperatorExpression::undefined;
00251 
00252   LinearListType::const_iterator i = __linearList.begin();
00253   while (i != __linearList.end()) {
00254     const IntegratedOperatorExpression& I = *(*i)->integrated();
00255 
00256     const IntegratedExpression& IE = *I.integratedExpression();
00257 
00258     if (IE.integratedExpressionType()
00259         == IntegratedExpression::function) {
00260       t = I.operatorType();
00261     }
00262     ++i;
00263   }
00264 
00265   return t;
00266 }
00267 
00268 VariationalOperator::Property
00269 MultiLinearExpression::
00270 getUnknownProperty() const
00271 {
00272   LinearListType::const_iterator i = __linearList.begin();
00273   while (i != __linearList.end()) {
00274     const IntegratedOperatorExpression& I = *(*i)->integrated();
00275 
00276     const IntegratedExpression& IE = *I.integratedExpression();
00277 
00278     if (IE.integratedExpressionType()
00279         == IntegratedExpression::unknownFunction) {
00280       if (I.isJump()) {
00281         return VariationalOperator::jump;
00282       }
00283       if (I.isMean()) {
00284         return VariationalOperator::mean;
00285       }
00286     }
00287     ++i;
00288   }
00289 
00290   return VariationalOperator::normal;
00291 }
00292 
00293 VariationalOperator::Property
00294 MultiLinearExpression::
00295 getTestFunctionProperty() const
00296 {
00297   LinearListType::const_iterator i = __linearList.begin();
00298   while (i != __linearList.end()) {
00299     const IntegratedOperatorExpression& I = *(*i)->integrated();
00300 
00301     const IntegratedExpression& IE = *I.integratedExpression();
00302 
00303     if (IE.integratedExpressionType()
00304         == IntegratedExpression::testFunction) {
00305       if (I.isJump()) {
00306         return VariationalOperator::jump;
00307       }
00308       if (I.isMean()) {
00309         return VariationalOperator::mean;
00310       }
00311     }
00312     ++i;
00313   }
00314 
00315   return VariationalOperator::normal;
00316 }
00317 
00318 ReferenceCounting<FunctionExpression>
00319 MultiLinearExpression::
00320 getFunctionExpression() const
00321 {
00322   ReferenceCounting<FunctionExpression> f = 0;
00323   LinearListType::const_iterator i = __linearList.begin();
00324   while (i != __linearList.end()) {
00325     const IntegratedExpression& IE =
00326       *(*i)->integrated()->integratedExpression();
00327 
00328     if (IE.integratedExpressionType()
00329         == IntegratedExpression::function) {
00330       const IntegratedExpressionFunctionExpression& I
00331         = dynamic_cast<const IntegratedExpressionFunctionExpression&>(IE);
00332 
00333       f = I.function();
00334     }
00335     ++i;
00336   }
00337   ASSERT(f != 0);
00338   return f;
00339 }
00340 
00341 MultiLinearExpression::OperatorType
00342 MultiLinearExpression::
00343 operatorType() const
00344 {
00345   IntegratedOperatorExpression::OperatorType unknown
00346     = IntegratedOperatorExpression::undefined;
00347 
00348   IntegratedOperatorExpression::OperatorType test
00349     = IntegratedOperatorExpression::undefined;
00350 
00351   IntegratedOperatorExpression::OperatorType function
00352     = IntegratedOperatorExpression::undefined;
00353 
00354   for(LinearListType::const_iterator i = __linearList.begin();
00355       i != __linearList.end(); ++i) {
00356     const IntegratedOperatorExpression& I = *(*i)->integrated();
00357     switch(I.integratedExpression()->integratedExpressionType()) {
00358     case IntegratedExpression::testFunction: {
00359       test = I.operatorType();
00360       break;
00361     }
00362     case IntegratedExpression::unknownFunction: {
00363       unknown = I.operatorType();
00364       break;
00365     }
00366     case IntegratedExpression::function: {
00367       function = I.operatorType();
00368       break;
00369     }
00370     }
00371   }
00372 
00373   OperatorType t = undefined;
00374 
00375   switch (__linearList.size()) {
00376   case 1: {
00377     switch (test) {
00378     case IntegratedOperatorExpression::orderZero: {
00379       t=V;
00380       break;
00381     }
00382     case IntegratedOperatorExpression::gradient: {
00383       t=gradV;
00384       break;
00385     }
00386     case IntegratedOperatorExpression::dx:
00387     case IntegratedOperatorExpression::dy:
00388     case IntegratedOperatorExpression::dz: {
00389       return dxV;
00390       break;
00391     }
00392     default: {
00393       throw ErrorHandler(__FILE__,__LINE__,
00394                          "not implemented",
00395                          ErrorHandler::unexpected);
00396     }
00397     }
00398     break;
00399   }
00400   case 2: {
00401     if (test != IntegratedOperatorExpression::undefined) {
00402       if (unknown != IntegratedOperatorExpression::undefined) {
00403         switch (unknown) {
00404         case IntegratedOperatorExpression::orderZero: {
00405           switch (test) {
00406           case IntegratedOperatorExpression::orderZero: {
00407             t=UV;
00408             break;
00409           }
00410           case IntegratedOperatorExpression::gradient:
00411           case IntegratedOperatorExpression::dx:
00412           case IntegratedOperatorExpression::dy:
00413           case IntegratedOperatorExpression::dz: {
00414             t=UdxV;
00415             break;
00416           }
00417           default: {
00418             throw ErrorHandler(__FILE__,__LINE__,
00419                                "not implemented",
00420                                ErrorHandler::unexpected);
00421           }
00422           }
00423           break;
00424         }
00425         case IntegratedOperatorExpression::gradient: {
00426           switch (test) {
00427           case IntegratedOperatorExpression::orderZero: {
00428             t=dxUV;
00429             break;
00430           }
00431           case IntegratedOperatorExpression::gradient: {
00432             t=gradUgradV;
00433             break;
00434           }
00435           case IntegratedOperatorExpression::dx:
00436           case IntegratedOperatorExpression::dy:
00437           case IntegratedOperatorExpression::dz: {
00438             return dxUdxV;
00439             break;
00440           }
00441           default: {
00442             throw ErrorHandler(__FILE__,__LINE__,
00443                                "not implemented",
00444                                ErrorHandler::unexpected);
00445           }
00446           }
00447           break;
00448         }
00449         case IntegratedOperatorExpression::dx:
00450         case IntegratedOperatorExpression::dy:
00451         case IntegratedOperatorExpression::dz: {
00452           switch (test) {
00453           case IntegratedOperatorExpression::orderZero: {
00454             t=dxUV;
00455             break;
00456           }
00457           case IntegratedOperatorExpression::gradient:
00458           case IntegratedOperatorExpression::dx:
00459           case IntegratedOperatorExpression::dy:
00460           case IntegratedOperatorExpression::dz: {
00461             t=dxUdxV;
00462             break;
00463           }
00464           default: {
00465             throw ErrorHandler(__FILE__,__LINE__,
00466                                "not implemented",
00467                                ErrorHandler::unexpected);
00468           }
00469           }
00470           break;
00471         }
00472         default: {
00473           throw ErrorHandler(__FILE__,__LINE__,
00474                              "not implemented",
00475                              ErrorHandler::unexpected);
00476         }
00477         }
00478       } else {
00479         switch (function) {
00480         case IntegratedOperatorExpression::gradient: {
00481           switch (test) {
00482           case IntegratedOperatorExpression::gradient: {
00483             t = gradFgradV;
00484             break;
00485           }
00486           case IntegratedOperatorExpression::orderZero:
00487           case IntegratedOperatorExpression::dx:
00488           case IntegratedOperatorExpression::dy:
00489           case IntegratedOperatorExpression::dz:
00490           default: {
00491             throw ErrorHandler(__FILE__,__LINE__,
00492                                "not implemented",
00493                                ErrorHandler::unexpected);
00494           }
00495           }
00496           break;
00497         }
00498         case IntegratedOperatorExpression::dx:
00499         case IntegratedOperatorExpression::dy:
00500         case IntegratedOperatorExpression::dz: {
00501           switch (test) {
00502           case IntegratedOperatorExpression::orderZero: {
00503             t = dxFV;
00504             break;
00505           }
00506           case IntegratedOperatorExpression::gradient:
00507           case IntegratedOperatorExpression::dx:
00508           case IntegratedOperatorExpression::dy:
00509           case IntegratedOperatorExpression::dz: {
00510             throw ErrorHandler(__FILE__,__LINE__,
00511                                "not implemented",
00512                                ErrorHandler::unexpected);
00513             break;
00514           }
00515           default: {
00516             throw ErrorHandler(__FILE__,__LINE__,
00517                                "not implemented",
00518                                ErrorHandler::unexpected);
00519           }
00520           }
00521           break;
00522         }
00523         default: {
00524           throw ErrorHandler(__FILE__,__LINE__,
00525                              "not implemented",
00526                              ErrorHandler::unexpected);
00527         }
00528         }
00529       }
00530     }
00531     break;
00532   }
00533   default: {
00534     throw ErrorHandler(__FILE__,__LINE__,
00535                        "not implemented",
00536                        ErrorHandler::unexpected);
00537   }
00538   }
00539   return t;
00540 }
00541 
00542 
00543 MultiLinearExpression::LinearFormType
00544 MultiLinearExpression::
00545 linearFormType() const
00546 {
00547   LinearFormType t;
00548   LinearListType::const_iterator i = __linearList.begin();
00549   switch (__linearList.size()) {
00550   case 1: {
00551     const IntegratedExpression& I=*(*i)->integrated()->integratedExpression();
00552     if (I.integratedExpressionType() != IntegratedExpression::testFunction) {
00553       const std::string errorMsg
00554         = "error defining "+stringify(*this)+"\n"
00555         "1-linear operator should be defined using test functions";
00556 
00557       throw ErrorHandler(__FILE__,__LINE__,
00558                          errorMsg,
00559                          ErrorHandler::normal);
00560     } else {
00561       return linear;
00562     }
00563     break;
00564   }
00565   case 2: {
00566     return biLinear;
00567     break;
00568   }
00569   default: {
00570     throw ErrorHandler(__FILE__,__LINE__,
00571                        "not implemented",
00572                        ErrorHandler::unexpected);
00573   }
00574   }
00575 
00576   return t;
00577 }
00578 
00579 
00580 void
00581 MultiLinearExpression::
00582 check()
00583 {
00584   if (__linearList.size() > 2) {
00585     this->execute();
00586     const std::string errorMsg
00587       = stringify(*this)+" is "+stringify(__linearList.size())
00588       +"-linear.\ncan discretize at most 2-linear expressions\n";
00589 
00590     throw ErrorHandler(__FILE__,__LINE__,
00591                        errorMsg,
00592                        ErrorHandler::normal);
00593   } else if (__linearList.size() == 2) {
00594     LinearListType::iterator i = __linearList.begin();
00595     IntegratedExpression::IType itype
00596       = (*(*(*(*i)).integrated()).integratedExpression()).integratedExpressionType();
00597     i++;
00598     IntegratedExpression::IType jtype
00599       = (*(*(*(*i)).integrated()).integratedExpression()).integratedExpressionType();
00600     if (itype == jtype) {
00601       this->execute();
00602       const std::string errorMsg 
00603         = stringify(*this)+" is non linear\n"
00604         +"Cannot discretize non linear operators\n";
00605 
00606       throw ErrorHandler(__FILE__,__LINE__,
00607                          errorMsg,
00608                          ErrorHandler::normal);
00609     }
00610   }
00611 }