"Fossies" - the Fresh Open Source Software Archive

Member "getdp-3.1.0-source/Functions/BF_Edge.cpp" (4 Apr 2019, 12588 Bytes) of package /linux/privat/getdp-3.1.0-source.tgz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "BF_Edge.cpp" see the Fossies "Dox" file reference documentation.

    1 // GetDP - Copyright (C) 1997-2019 P. Dular and C. Geuzaine, University of Liege
    2 //
    3 // See the LICENSE.txt file for license information. Please report all
    4 // issues on https://gitlab.onelab.info/getdp/getdp/issues.
    5 //
    6 // Contributor(s):
    7 //   Christophe Trophime
    8 //
    9 
   10 #include "ProData.h"
   11 #include "Message.h"
   12 
   13 #define SQU(a)     ((a)*(a))
   14 
   15 #define NoEdge  Message::Error("Missing Edge Entity in Element %d", Element->Num)
   16 
   17 /* ------------------------------------------------------------------------ */
   18 /*  B F _ E d g e                                                           */
   19 /* ------------------------------------------------------------------------ */
   20 
   21 #define WrongNumEdge   Message::Error("Wrong Edge number in 'BF_Edge'")
   22 
   23 void BF_Edge(struct Element * Element, int NumEdge,
   24          double u, double v, double w,  double s[])
   25 {
   26   switch (Element->Type) {
   27   case LINE :
   28     switch(NumEdge) {
   29     case 1  : s[0] = 0.5 ; s[1] = 0. ; s[2] = 0. ; break ;
   30     default : WrongNumEdge ;
   31     }
   32     break ;
   33 
   34   case TRIANGLE :
   35   case TRIANGLE_2 :
   36     switch(NumEdge) {
   37     case 1  : s[0] = 1.-v ; s[1] = u    ; s[2] = 0.  ; break ;
   38     case 2  : s[0] = v    ; s[1] = 1.-u ; s[2] = 0.  ; break ;
   39     case 3  : s[0] =   -v ; s[1] =    u ; s[2] = 0.  ; break ;
   40     default : WrongNumEdge ;
   41     }
   42     break ;
   43 
   44   case QUADRANGLE :
   45     switch(NumEdge) {
   46     case 1  : s[0] =  0.25 * (1.-v) ; s[1] = 0.            ; s[2] = 0. ; break ;
   47     case 2  : s[0] =  0.            ; s[1] = 0.25 * (1.-u) ; s[2] = 0. ; break ;
   48     case 3  : s[0] =  0.            ; s[1] = 0.25 * (1.+u) ; s[2] = 0. ; break ;
   49     case 4  : s[0] = -0.25 * (1.+v) ; s[1] = 0.            ; s[2] = 0. ; break ;
   50     default : WrongNumEdge ;
   51     }
   52     break ;
   53 
   54   case TETRAHEDRON :
   55     switch(NumEdge) {
   56     case 1  : s[0] =  1.-v-w ; s[1] =  u      ; s[2] = u      ; break ;
   57     case 2  : s[0] =  v      ; s[1] =  1.-u-w ; s[2] = v      ; break ;
   58     case 3  : s[0] =  w      ; s[1] =  w      ; s[2] = 1.-u-v ; break ;
   59     case 4  : s[0] = -v      ; s[1] =  u      ; s[2] = 0.     ; break ;
   60     case 5  : s[0] = -w      ; s[1] =  0.     ; s[2] = u      ; break ;
   61     case 6  : s[0] =  0.     ; s[1] = -w      ; s[2] = v      ; break ;
   62     default : WrongNumEdge ;
   63     }
   64     break ;
   65 
   66   case HEXAHEDRON :
   67     switch(NumEdge) {
   68     case 1  : s[0] =  0.125 * (1.-v) * (1.-w) ; s[1] = 0. ; s[2] = 0. ; break ;
   69     case 6  : s[0] = -0.125 * (1.+v) * (1.-w) ; s[1] = 0. ; s[2] = 0. ; break ;
   70     case 9  : s[0] =  0.125 * (1.-v) * (1.+w) ; s[1] = 0. ; s[2] = 0. ; break ;
   71     case 12 : s[0] = -0.125 * (1.+v) * (1.+w) ; s[1] = 0. ; s[2] = 0. ; break ;
   72 
   73     case 2  : s[0] = 0. ; s[1] = 0.125 * (1.-u) * (1.-w)  ; s[2] = 0. ; break ;
   74     case 4  : s[0] = 0. ; s[1] = 0.125 * (1.+u) * (1.-w)  ; s[2] = 0. ; break ;
   75     case 10 : s[0] = 0. ; s[1] = 0.125 * (1.-u) * (1.+w)  ; s[2] = 0. ; break ;
   76     case 11 : s[0] = 0. ; s[1] = 0.125 * (1.+u) * (1.+w)  ; s[2] = 0. ; break ;
   77 
   78     case 3  : s[0] = 0. ; s[1] = 0. ; s[2] = 0.125 * (1.-u) * (1.-v)  ; break ;
   79     case 5  : s[0] = 0. ; s[1] = 0. ; s[2] = 0.125 * (1.+u) * (1.-v)  ; break ;
   80     case 7  : s[0] = 0. ; s[1] = 0. ; s[2] = 0.125 * (1.+u) * (1.+v)  ; break ;
   81     case 8  : s[0] = 0. ; s[1] = 0. ; s[2] = 0.125 * (1.-u) * (1.+v)  ; break ;
   82     default : WrongNumEdge ;
   83     }
   84     break ;
   85 
   86   case PRISM :
   87     switch(NumEdge) {
   88     case 1  : s[0] =  0.5 * (1.-v) * (1.-w) ;
   89               s[1] =  0.5 * u      * (1.-w) ;
   90               s[2] =  0.                    ; break ;
   91     case 2  : s[0] =  0.5 * v      * (1.-w) ;
   92               s[1] =  0.5 * (1.-u) * (1.-w) ;
   93               s[2] =  0.                    ; break ;
   94     case 3  : s[0] =  0.                    ;
   95               s[1] =  0.                    ;
   96               s[2] =  0.5 * (1.-u-v)        ; break ;
   97     case 4  : s[0] = -0.5 * v      * (1.-w) ;
   98               s[1] =  0.5 * u      * (1.-w) ;
   99               s[2] =  0.                    ; break ;
  100     case 5  : s[0] =  0.                    ;
  101               s[1] =  0.                    ;
  102               s[2] =  0.5 * u               ; break ;
  103     case 6  : s[0] =  0.                    ;
  104               s[1] =  0.                    ;
  105               s[2] =  0.5 * v               ; break ;
  106     case 7  : s[0] =  0.5 * (1.-v) * (1.+w) ;
  107               s[1] =  0.5 * u      * (1.+w) ;
  108               s[2] =  0.                    ; break ;
  109     case 8  : s[0] =  0.5 * v      * (1.+w) ;
  110               s[1] =  0.5 * (1.-u) * (1.+w) ;
  111               s[2] =  0.                    ; break ;
  112     case 9  : s[0] = -0.5 * v      * (1.+w) ;
  113               s[1] =  0.5 * u      * (1.+w) ;
  114               s[2] =  0.                    ; break ;
  115     default : WrongNumEdge ;
  116     }
  117     break ;
  118 
  119   case PYRAMID :
  120     if (w != 1){
  121       switch(NumEdge) {
  122       case 1  : s[0] =  0.25 * (1 - v - w) ;
  123             s[1] =  0. ;
  124         s[2] =  0.25 * (u - u * v / (1. - w)) ; break ;
  125       case 2  : s[0] =  0. ;
  126             s[1] =  0.25 * (1 - u - w) ;
  127         s[2] =  0.25 * (v - u * v / (1. - w))  ; break ;
  128       case 4  : s[0] =  0. ;
  129                 s[1] =  0.25 * (1 + u - w) ;
  130         s[2] =  0.25 * (v + u * v / (1. - w)) ; break ;
  131       case 6  : s[0] = -0.25 * (1 + v - w) ;
  132                 s[1] =  0. ;
  133         s[2] = -0.25 * (u + u * v / (1. - w)) ; break ;
  134       case 3  : s[0] =  0.25 * (w - v * w / (1. - w)) ;
  135                 s[1] =  0.25 * (w - u * w / (1. - w)) ;
  136         s[2] =  0.25 * (1. - u - v + u * v / SQU(1. - w) - 2 * u * v * w / SQU(1. - w)) ; break ;
  137       case 5  : s[0] = -0.25 * (w - v * w / (1. - w)) ;
  138                 s[1] =  0.25 * (w + u * w / (1. - w)) ;
  139         s[2] =  0.25 * (1. + u - v - u * v / SQU(1. - w) + 2 * u * v * w / SQU(1. - w)) ; break ;
  140       case 7  : s[0] = -0.25 * (w + v * w / (1. - w)) ;
  141                 s[1] = -0.25 * (w + u * w / (1. - w)) ;
  142                 s[2] =  0.25 * (1. + u + v + u * v / SQU(1. - w) - 2 * u * v * w / SQU(1. - w)) ; break ;
  143       case 8  : s[0] =  0.25 * (w + v * w / (1. - w)) ;
  144             s[1] = -0.25 * (w - u * w / (1. - w)) ;
  145         s[2] =  0.25 * (1. - u + v - u * v / SQU(1. - w) + 2 * u * v * w / SQU(1. - w)) ; break ;
  146       default : WrongNumEdge ;
  147       }
  148     }
  149     else
  150       switch(NumEdge) {
  151       case 1  : s[0] = -0.25 * v ;
  152             s[1] =  0. ;
  153         s[2] =  0.25 * u ; break ;
  154       case 2  : s[0] =  0. ;
  155             s[1] = -0.25 * u ;
  156         s[2] =  0.25 * v ; break ;
  157       case 4  : s[0] =  0. ;
  158                 s[1] =  0.25 * u ;
  159         s[2] =  0.25 * v ; break ;
  160       case 6  : s[0] = -0.25 * v ;
  161                 s[1] =  0. ;
  162         s[2] = -0.25 * u ; break ;
  163       case 3  : s[0] =  0.25 ;
  164             s[1] =  0.25 ;
  165         s[2] =  0.25 * (1. - u - v) ; break ;
  166       case 5  : s[0] = -0.25 ;
  167                 s[1] =  0.25 ;
  168         s[2] =  0.25 * (1. + u - v) ; break ;
  169       case 7  : s[0] = -0.25 ;
  170                 s[1] = -0.25 ;
  171                 s[2] =  0.25 * (1. + u + v) ; break ;
  172       case 8  : s[0] =  0.25 ;
  173             s[1] = -0.25 ;
  174         s[2] =  0.25 * (1. - u + v) ; break ;
  175       default : WrongNumEdge ;
  176       }
  177     break ;
  178 
  179   default :
  180     Message::Error("Unknown type of Element in BF_Edge");
  181     break;
  182   }
  183 
  184   if (!Element->GeoElement->NumEdges) NoEdge ;
  185 
  186   if (Element->GeoElement->NumEdges[NumEdge-1] < 0) {
  187     s[0] = - s[0] ; s[1] = - s[1] ; s[2] = - s[2] ;
  188   }
  189 }
  190 
  191 #undef WrongNumEdge
  192 
  193 /* ------------------------------------------------------------------------ */
  194 /*  B F _ C u r l E d g e                                                   */
  195 /* ------------------------------------------------------------------------ */
  196 
  197 #define WrongNumEdge   Message::Error("Wrong Edge number in 'BF_CurlEdge'")
  198 
  199 void BF_CurlEdge(struct Element * Element, int NumEdge,
  200          double u, double v, double w,  double s[])
  201 {
  202   switch (Element->Type) {
  203   case LINE :
  204     switch(NumEdge) {
  205     case 1  : s[0] = 0. ; s[1] = 0. ; s[2] = 0. ; break ;
  206     default : WrongNumEdge ;
  207     }
  208     break ;
  209 
  210   case TRIANGLE :
  211   case TRIANGLE_2 :
  212     switch(NumEdge) {
  213     case 1  : s[0] = 0. ; s[1] = 0. ; s[2] =  2. ; break ;
  214     case 2  : s[0] = 0. ; s[1] = 0. ; s[2] = -2. ; break ;
  215     case 3  : s[0] = 0. ; s[1] = 0. ; s[2] =  2. ; break ;
  216     default : WrongNumEdge ;
  217     }
  218     break ;
  219 
  220   case QUADRANGLE :
  221     switch(NumEdge) {
  222     case 1  : s[0] = 0. ; s[1] = 0. ; s[2] =  0.25 ; break ;
  223     case 2  : s[0] = 0. ; s[1] = 0. ; s[2] = -0.25 ; break ;
  224     case 3  : s[0] = 0. ; s[1] = 0. ; s[2] =  0.25 ; break ;
  225     case 4  : s[0] = 0. ; s[1] = 0. ; s[2] =  0.25 ; break ;
  226     default : WrongNumEdge ;
  227     }
  228     break ;
  229 
  230   case TETRAHEDRON :
  231     switch(NumEdge) {
  232     case 1  : s[0] =  0. ; s[1] = -2. ; s[2] =  2. ; break ;
  233     case 2  : s[0] =  2. ; s[1] =  0. ; s[2] = -2. ; break ;
  234     case 3  : s[0] = -2. ; s[1] =  2. ; s[2] =  0. ; break ;
  235     case 4  : s[0] =  0. ; s[1] =  0. ; s[2] =  2. ; break ;
  236     case 5  : s[0] =  0. ; s[1] = -2. ; s[2] =  0. ; break ;
  237     case 6  : s[0] =  2. ; s[1] =  0. ; s[2] =  0. ; break ;
  238     default : WrongNumEdge ;
  239     }
  240     break ;
  241 
  242   case HEXAHEDRON :
  243     switch(NumEdge) {
  244     case 1  : s[0] = 0. ; s[1] = 0.125*(v-1.) ; s[2] = 0.125*(1.-w) ; break ;
  245     case 6  : s[0] = 0. ; s[1] = 0.125*(v+1.) ; s[2] = 0.125*(1.-w) ; break ;
  246     case 9  : s[0] = 0. ; s[1] = 0.125*(1.-v) ; s[2] = 0.125*(1.+w) ; break ;
  247     case 12 : s[0] = 0. ; s[1] =-0.125*(v+1.) ; s[2] = 0.125*(1.+w) ; break ;
  248 
  249     case 2  : s[0] = 0.125*(1.-u) ; s[1] = 0. ; s[2] = 0.125*(w-1.) ; break ;
  250     case 4  : s[0] = 0.125*(1.+u) ; s[1] = 0. ; s[2] = 0.125*(1.-w) ; break ;
  251     case 10 : s[0] = 0.125*(u-1.) ; s[1] = 0. ; s[2] =-0.125*(w+1.) ; break ;
  252     case 11 : s[0] =-0.125*(1.+u) ; s[1] = 0. ; s[2] = 0.125*(w+1.) ; break ;
  253 
  254     case 3  : s[0] = 0.125*(u-1.) ; s[1] = 0.125*(1.-v) ; s[2] = 0. ; break ;
  255     case 5  : s[0] =-0.125*(u+1.) ; s[1] = 0.125*(v-1.) ; s[2] = 0. ; break ;
  256     case 7  : s[0] = 0.125*(u+1.) ; s[1] =-0.125*(1.+v) ; s[2] = 0. ; break ;
  257     case 8  : s[0] = 0.125*(1.-u) ; s[1] = 0.125*(1.+v) ; s[2] = 0. ; break ;
  258     default : WrongNumEdge ;
  259     }
  260     break ;
  261 
  262   case PRISM :
  263     switch(NumEdge) {
  264     case 1  : s[0] =  0.5*u      ; s[1] =  0.5*(v-1.) ; s[2] =  1.-w ; break ;
  265     case 2  : s[0] =  0.5*(1.-u) ; s[1] = -0.5*v      ; s[2] =  w-1. ; break ;
  266     case 3  : s[0] = -0.5        ; s[1] =  0.5        ; s[2] =  0.   ; break ;
  267 
  268     case 4  : s[0] =  0.5*u      ; s[1] =  0.5*v      ; s[2] =  1.-w ; break ;
  269     case 5  : s[0] =  0.         ; s[1] = -0.5        ; s[2] =  0.   ; break ;
  270     case 6  : s[0] =  0.5        ; s[1] =  0.         ; s[2] =  0.   ; break ;
  271 
  272     case 7  : s[0] = -0.5*u      ; s[1] =  0.5*(1.-v) ; s[2] =  1.+w ; break ;
  273     case 8  : s[0] =  0.5*(u-1.) ; s[1] =  0.5*v      ; s[2] = -1.-w ; break ;
  274     case 9  : s[0] = -0.5*u      ; s[1] = -0.5*v      ; s[2] =  1.+w ; break ;
  275     default : WrongNumEdge ;
  276     }
  277     break ;
  278 
  279   case PYRAMID :
  280     if (w != 1){
  281       switch(NumEdge) {
  282       case 1  : s[0] = -0.25 * u / (1. - w) ;       s[1] = -0.5 + 0.25 * v / (1. - w) ; s[2] =  0.25 ; break ;
  283       case 2  : s[0] =  0.5 - 0.25 * u / (1. - w) ; s[1] =  0.25 * v / (1. - w) ;       s[2] = -0.25 ; break ;
  284       case 4  : s[0] =  0.5 + 0.25 * u / (1. - w) ; s[1] = -0.25 * v / (1. - w) ;       s[2] =  0.25 ; break ;
  285       case 6  : s[0] = -0.25 * u / (1. - w) ;       s[1] =  0.5 + 0.25 * v / (1. - w) ; s[2] =  0.25 ; break ;
  286       case 3  : s[0] = -0.5 * (1. - u / (1. - w)) ; s[1] =  0.5 * (1. - v / (1. - w)) ; s[2] =  0. ; break;
  287       case 5  : s[0] = -0.5 * (1. + u / (1. - w)) ; s[1] = -0.5 * (1. - v / (1. - w)) ; s[2] =  0. ; break;
  288       case 7  : s[0] =  0.5 * (1. + u / (1. - w)) ; s[1] = -0.5 * (1. + v / (1. - w)) ; s[2] =  0. ; break;
  289       case 8  : s[0] =  0.5 * (1. - u / (1. - w)) ; s[1] =  0.5 * (1. + v / (1. - w)) ; s[2] =  0. ; break;
  290       default : WrongNumEdge ;
  291       }
  292     } else {
  293       switch(NumEdge) {
  294       case 1  : s[0] =  0.  ; s[1] = -0.5 ; s[2] =  0.25 ; break ;
  295       case 2  : s[0] =  0.5 ; s[1] =  0.  ; s[2] = -0.25 ; break ;
  296       case 4  : s[0] =  0.5 ; s[1] =  0.  ; s[2] =  0.25 ; break ;
  297       case 6  : s[0] =  0.  ; s[1] =  0.5 ; s[2] =  0.25 ; break ;
  298       case 3  : s[0] = -0.5 ; s[1] =  0.5 ; s[2] =  0. ; break;
  299       case 5  : s[0] = -0.5 ; s[1] = -0.5 ; s[2] =  0. ; break;
  300       case 7  : s[0] =  0.5 ; s[1] = -0.5 ; s[2] =  0. ; break;
  301       case 8  : s[0] =  0.5 ; s[1] =  0.5 ; s[2] =  0. ; break;
  302       default : WrongNumEdge ;
  303       }
  304     }
  305     break ;
  306 
  307 
  308   default :
  309     Message::Error("Unknown type of Element in BF_CurlEdge");
  310     break;
  311   }
  312 
  313   if (!Element->GeoElement->NumEdges) NoEdge ;
  314 
  315   if (Element->GeoElement->NumEdges[NumEdge-1] < 0) {
  316     s[0] = - s[0] ; s[1] = - s[1] ; s[2] = - s[2] ;
  317   }
  318 }
  319 
  320 #undef WrongNumEdge
  321 
  322 #undef NoEdge