domain_bispheric_rect_import.cpp

/*
    Copyright 2017 Philippe Grandclement
 
    This file is part of Kadath.
 
    Kadath is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
 
    Kadath is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.
 
    You should have received a copy of the GNU General Public License
    along with Kadath.  If not, see <http://www.gnu.org/licenses/>.
*/
 
#include "headcpp.hpp"
 
#include "utilities.hpp"
#include "bispheric.hpp"
#include "term_eq.hpp"
#include "scalar.hpp"
#include "tensor_impl.hpp"
 
namespace Kadath {
// Tensorial parts :
Tensor Domain_bispheric_rect::import (int numdom, int bound, int n_ope, const Array<int>& zedoms, Tensor** parts) const {
  
  Tensor res (*parts[0], false);
    
  for (int nc=0 ; nc<res.get_n_comp() ; nc++)
    res.cmp[nc]->set_domain(numdom).allocate_conf() ;
  
   // Need coefficients of parts :
  for (int i=0 ; i<n_ope ; i++)
    for (int nc=0 ; nc<parts[i]->get_n_comp() ; nc++)
    parts[i]->cmp[nc]->set_domain(zedoms(i)).coef() ;
 
  // Loop on the points of the boundary :
  Val_domain xx (get_cart(1)) ;
  Val_domain yy (get_cart(2)) ;
  Val_domain zz (get_cart(3)) ;
  
  if (bound==INNER_BC) {
      int index_mu = 0 ;
      Index pos (get_nbr_points()) ;
      Index pos_bound (get_nbr_points()) ;
  
      for (int k=0 ; k<nbr_points(2) ; k++)
    for (int j=0 ; j<nbr_points(1) ; j++) {
     
    // Indices in the shell
    pos_bound.set(0) = index_mu ;
    pos_bound.set(1) = j ;
    pos_bound.set(2) = k ;
      
    // Absolute coordinates
    Point MM (3) ;
    MM.set(1) = xx(pos_bound) ;
    MM.set(2) = yy(pos_bound) ;
    MM.set(3) = zz(pos_bound) ;
      
    // In which other domain is it ?
    bool found = false ;
    int current = 0 ;
    while ((current<n_ope) && (!found)) {
      if (parts[0]->get_space().get_domain(zedoms(current))->is_in(MM))
        found = true ;
    else 
      current ++ ;
      }
      if (!found) {
    cerr << "Point " << MM << " not found in other domains, for Domain_bispheric_rect::import" << endl ;
    abort() ;
      }
      
      // Convert to numerical coordinates of the other domain
      Point num (parts[0]->get_space().get_domain(zedoms(current))->absol_to_num(MM)) ;
      
      // Now loop on the components :
      for (int nc=0 ; nc<res.get_n_comp() ; nc++) {
    double val = ((*parts[current]->cmp[nc])(zedoms(current)).check_if_zero()) ? 0 : (*parts[current]->cmp[nc])(zedoms(current)).get_base().summation (num, *(*parts[current]->cmp[nc])(zedoms(current)).cf) ;
    // Loop on radius :
    for (int i=0 ; i<nbr_points(0) ; i++) {
      pos.set(0) = i ;
      pos.set(1) = j ;
      pos.set(2) = k ;
     
      res.cmp[nc]->set_domain(numdom).set(pos) = val ;
    }
      }
    }
    
    // Assert a std_base :
    res.std_base() ;
  return res ;
  }
  
  
  if (bound==OUTER_BC) {
      int index_mu = nbr_points(0)-1 ;
      int index_chi = nbr_points(1)-1 ;
      
      Index pos (get_nbr_points()) ;
      Index pos_bound (get_nbr_points()) ;
  
      for (int k=0 ; k<nbr_points(2) ; k++) {
     
    // Indices in the shell
    pos_bound.set(0) = index_mu ;
    pos_bound.set(1) = index_chi ;
    pos_bound.set(2) = k ;
      
    // Absolute coordinates
    Point MM (3) ;
    MM.set(1) = xx(pos_bound) ;
    MM.set(2) = yy(pos_bound) ;
    MM.set(3) = zz(pos_bound) ;
      
    // In which other domain is it ?
    bool found = false ;
    int current = 0 ;
    while ((current<n_ope) && (!found)) {
      if (parts[0]->get_space().get_domain(zedoms(current))->is_in(MM))
        found = true ;
    else 
      current ++ ;
      }
      if (!found) {
    cerr << "Point " << MM << " not found in other domains, for Domain_bispheric_rect::import" << endl ;
    abort() ;
      }
      
      // Convert to numerical coordinates of the other domain
      Point num (parts[0]->get_space().get_domain(zedoms(current))->absol_to_num(MM)) ;
      
      // Now loop on the components :
      for (int nc=0 ; nc<res.get_n_comp() ; nc++) {
    double val = ((*parts[current]->cmp[nc])(zedoms(current)).check_if_zero()) ? 0 : (*parts[current]->cmp[nc])(zedoms(current)).get_base().summation (num, *(*parts[current]->cmp[nc])(zedoms(current)).cf) ;
    // Loop on radius :
    for (int i=0 ; i<nbr_points(0) ; i++) 
      for (int j=0 ; j<nbr_points(1) ; j++) {
      pos.set(0) = i ;
      pos.set(1) = j ;
      pos.set(2) = k ;
     
      res.cmp[nc]->set_domain(numdom).set(pos) = val ;
    }
      }
    }
        
    // Assert a std_base :
    res.std_base() ;
    return res ;
  }
  
  cerr << "unknown bound in Domain_bispheric_rect::import" << endl ;
  abort() ;    
}}