meer.cpp

Zobrazit dokumentaci tohoto souboru.

#include "meer.h"

#ifdef DEVEL
#include "arrays.h"
#include "gpa.h"
#else
#include "arrays.h"
#include "gpa.h"
#endif

#include <stdio.h>
#include <stdlib.h>
#include <string.h>


namespace meerspace {

MEER :: MEER (void)
{
  // private
  refined_node_values_1 = refined_node_values_2 = NULL;
  valtype_1 = valtype_2 = MEER_IPV_Void;
  
  // protected
  enorm = MEER_EN_Void;
  h_max_ratio = 2.0;    // max zvetsit na X-nasobek
  h_min_ratio = 0.25;   // max zmensit na X-nasobek
  SPR_boundary = MEER_SPRbpt_basic;
  
  normal_at_nodes = false;
  
  nREGIONs = 0;  REGIONs = NULL;
  nNODEs   = 0;  NODEs = NULL;
  nELEMs   = 0;  ELEMs = NULL;
}

MEER :: ~MEER ()
{
  // private
  long i, j;
  if (refined_node_values_1) {
    for (i=0; i<nREGIONs; i++) {
      for (j=0; j<nNODEs; j++)
    delete refined_node_values_1[i][j];
      
      delete [] refined_node_values_1[i];
    }
    delete [] refined_node_values_1;
  }
  
  if (refined_node_values_2) {
    for (i=0; i<nREGIONs; i++) {
      for (j=0; j<nNODEs; j++)
    delete refined_node_values_2[i][j];
      
      delete [] refined_node_values_2[i];
    }
    delete [] refined_node_values_2;
  }
  
  // protected
  for (i=0; i<nNODEs; i++) {
    if (NODEs[i].normal == NULL)
      errol;
    for (j=0; j<nREGIONs; j++)
      delete NODEs[i].normal[j];
  }
  
  delete [] REGIONs;
  delete [] NODEs;
  delete [] ELEMs;
}


void MEER :: set_alloc_REGIONs (long n) { nREGIONs = n;  REGIONs = new REGION [nREGIONs]; for (long i=0; i<nREGIONs; i++)  REGIONs[i].id = i; } 
void MEER :: set_alloc_NODEs   (long n) { nNODEs   = n;  NODEs   = new NODE   [nNODEs];   for (long i=0; i<nNODEs; i++)    NODEs[i].id = i;   }
void MEER :: set_alloc_ELEMs   (long n) { nELEMs   = n;  ELEMs   = new ELEMENT[nELEMs];   for (long i=0; i<nELEMs; i++)    ELEMs[i].id = i;   }


//* main function
double MEER :: solve (double rerror)
{
  // boundary nodes
  this->set_boundary_nodes();
  
  //
  switch (enorm) {
  case MEER_EN_strain_stress:  valtype_1 = MEER_IPV_strain;  valtype_2 = MEER_IPV_stress;  break;
  case MEER_EN_stress_stress:  valtype_1 = MEER_IPV_stress;  valtype_2 = MEER_IPV_Void;    break;
  default: _errorr ("Unknown enorm");
  }
  
  // nodal recovery of values - z IP interpoluje do uzlu slozky zvolene veliciny => vysledkem je vylepsena/recovered velicina v uzlech
  this->MEER_SPR_nodal_recovery();
  
  // error computation
  this->MEER_error_estimatior ();
  
  // mesh refinement
  return this->MEER_mesh_refinement (rerror);
}

void MEER :: set_boundary_nodes (void)
{
  ;
}

void MEER :: MEER_error_estimatior (void)
{
  long j, k, l, regid, nip, nvalcomps, nno;
  long *nodes;
  double w, dV, V;
  long elem;
  Dvctr n(4), IPvalsR_1, IPvalsR_2, IPvalsF_1, IPvalsF_2;
  GPA<const Dvctr> recovered_nodal_values_1;
  GPA<const Dvctr> recovered_nodal_values_2;
  
  e2i.resize_ignore_vals(nELEMs);  e2i.zero();
  u2i.resize_ignore_vals(nELEMs);  u2i.zero();
  
  for (elem=0; elem<nELEMs; elem++) {
    nno       = ELEMs[elem].nnodes;
    nodes     = ELEMs[elem].nodes;
    regid     = ELEMs[elem].regid;
    nvalcomps = REGIONs[regid].nvals;
    
    nip = this->give_number_of_IPs2(elem);
    
    if (nno == 3 && nip < 3) _errorr("low number of int points");
    if (nno == 4 && nip < 4) _errorr("low number of int points");
    
    // assemble recovered/refined values at element nodes
    if (valtype_1) {
      recovered_nodal_values_1.resize_zero();
      for (j=0; j<nno; j++)
    recovered_nodal_values_1.add(refined_node_values_1[regid][nodes[j]]);
    }
    if (valtype_2) {
      recovered_nodal_values_2.resize_zero();
      for (j=0; j<nno; j++)
    recovered_nodal_values_2.add(refined_node_values_2[regid][nodes[j]]);
    }
    
    IPvalsR_1.resize_ignore_vals(nvalcomps);
    IPvalsR_2.resize_ignore_vals(nvalcomps);
    
    // element volume reduced
    V = this->give_element_area(elem) * this->give_element_thickness_reduced(elem);
    
    // error estimation at element - numericky, tj. smyce pres IPs
    for (j=0; j<nip; j++) {
      // IPvalsR - Refined values - ziskaji se interpolaci z uzlu do IP
      w = this->eval_interpol_fces_at_IPs2 (elem, j, n.give_ptr2val());
      
      if (valtype_1) {
    IPvalsR_1.zero();
    for (k=0; k<nvalcomps; k++)
      for (l=0; l<nno; l++)
        IPvalsR_1[k] += n[l] * recovered_nodal_values_1[l][0][k];
      }
      if (valtype_2) {
    IPvalsR_2.zero();
    for (k=0; k<nvalcomps; k++)
      for (l=0; l<nno; l++)
        IPvalsR_2[k] += n[l] * recovered_nodal_values_2[l][0][k];
      }
      
      // IPvalsF - FEM values - ziskane z FEM vypoctu v IP
      if (valtype_1)  IPvalsF_1.assign_array(this->give_IPs2_values(elem, j, valtype_1), nvalcomps);
      if (valtype_2)  IPvalsF_2.assign_array(this->give_IPs2_values(elem, j, valtype_2), nvalcomps);
      
      // error = IPvalsR - IPvalsF
      if (valtype_1)  IPvalsR_1.sbt(IPvalsF_1);
      if (valtype_2)  IPvalsR_2.sbt(IPvalsF_2);
      
      dV = w * V;
      //dV = w * this->give_element_volume(elem);
      
      switch (enorm) {
      case MEER_EN_strain_stress:
    e2i[elem] += dV * IPvalsR_1.give_dotproduct(IPvalsR_2);
    u2i[elem] += dV * IPvalsF_1.give_dotproduct(IPvalsF_2);
    break;
      case MEER_EN_stress_stress:
    e2i[elem] += dV * IPvalsR_1.give_dotproduct(IPvalsR_1);
    u2i[elem] += dV * IPvalsF_1.give_dotproduct(IPvalsF_1);
    break;
      default: _errorr ("Unsupported enorm");
      }
      
      if (valtype_1)  IPvalsF_1.free();
      if (valtype_2)  IPvalsF_2.free();
    }
  }
}

double MEER :: MEER_mesh_refinement (double rerror)
{
  // zaporna energie
  if (1) {
    long i,j,k,n;
    long nsuperelems, superelems[64];
    double sum;
    
    for (i=0; i<nELEMs; i++) {
      if (e2i[i] >= 0.0) continue;
      sum = 0.0;
      n = 0;
      for (j=0; j<ELEMs[i].nnodes; j++) {
    this->give_superelems_to_node(NODEs[ELEMs[i].nodes[j]].id, nsuperelems, superelems);
    for (k=0; k<nsuperelems; k++) {
      if (e2i[superelems[k]] >= 0.0) {
        sum += e2i[superelems[k]];
        n++;
      }
    }
      }
      e2i[i] = sum / n;
    }
  }
  
  double e2 = e2i.give_sum();
  double u2 = u2i.give_sum();
  
  double pe = sqrt( e2 / ( e2 + u2 ) );
  
  //if (true || pe > rerror) {
  long i;
  double uaver, em, coeff, ksi;
  
  // em = limit error at element
  coeff = 1.0;  // nevim proc to tady borek ma, ma tam 2.0
  em = coeff * rerror * sqrt( ( u2 + e2 ) / nELEMs );
  uaver =               sqrt( ( u2 + e2 ) / nELEMs );
  
  //hnewi_node.resize_ignore_vals(nNODEs);  hnewi_node.zero();
  h_old.resize_ignore_vals(nELEMs);  h_old.zero();
  h_new.resize_ignore_vals(nELEMs);  h_new.zero();
  
  errori_loc.resize_ignore_vals(nELEMs);
  errori_glob.resize_ignore_vals(nELEMs);
  
  for (i=0; i<nELEMs; i++) {
    // percentual error per element
    errori_loc [i] = sqrt( e2i[i]   / ( u2i[i] + e2i[i] ) ) * 100;
    errori_glob[i] = sqrt( e2i[i] ) /   uaver               * 100;
    
    ksi = sqrt(e2i[i]) / em;
    
    //
    if ( ksi < 1.0/h_max_ratio )  ksi = 1.0/h_max_ratio;  // prvek se zvetsi maximalne h_max krat
    if ( ksi > 1.0/h_min_ratio )  ksi = 1.0/h_min_ratio;  // prvek se zmensi maximalne h_min krat
    
    h_old[i] = this->give_element_characteristic_size(i);   // stara char. velikost prvku
    // tady bych mel zkontrolovat stupen aproximace, pokud to neni 1, tak se vzorecek zmeni na hnew = hold / pow(ksi, 1.0 / ord);
    h_new[i] = h_old[i] / ksi;
  }
  
  // zatim docasne, pro risovo hezky obrazek
  ei_aver.resize_ignore_vals(nELEMs);
  for (i=0; i<nELEMs; i++) {
    ei_aver[i] = e2i[i] / this->give_element_area(i) / this->give_element_thickness(i);
  }
  
  return pe;
}

// ******************************************************************************************
// - SPR_nodal_recovery == metoda, ktera z IP interpoluje do uzlu slozky zvolene veliciny
//   => vysledkem je vylepsena/recovered velicina v uzlech
// - do uzlu se interpoluje pomoci patche=zaplaty tvorene nekolika elementu k uzle prilehajich
//
// bp_appli == basic patch applicability
// =====================================
// - basic patch je zaplata z elementu, ktere uzel obsahuji (a jsou v danem regionu)
// - bp_appli nabyva hodnot:
//   0 - out of region          == patch je nulovy     - vylepsena velicina v uzlu se NEPOCITA 
//   1 - directly recovered     == patch je dostatecny - vylepsena velicina v uzlu se pocita ze zakladni patche
//   2 - in-directly recovered  == patch je maly       - vylepsena velicina v uzlu se pocita nejak jinak, treba z patche sousedniho uzlu s bpa==1
//   x - uz vic cisel nepridavej, pocita se s tim v kodu
// - bp_elems - pole ID elementu patrici do bp
//
// A) predpokladam, ze kazdy 2 lezi aspon v jedne patchi 1, pak se pri vypoctu 1 vycisli hodnota i na 2 [AKTUALNI RESENI]
//    2 muze patrit do vice patchi 1
//    a) muzu je vzit vsechny a zprumerovat [AKTUALNI RESENI}
//    b) muzu vzit 1, v ktere lezi vsechny superelems pro 2
//       - toto je asi presnejsi, ale ja na to zatim prdim, celkova nepresnost vypoctu je tak velka ze se mi s tim nechce delat
//         ale kdyz bude cas tak to udelam
// B) pro kazdy 2 vezmu superelemns plus jeste elementy, ktere se hranou dotykaji superelems
//    to by melo stacit pro dostatecnej pocet nsp
// C) oofem ma jeste variantu metody A), kdy vsechny boundary uzly jsou 2 a stavaji se 1 jen kdyz se znich ma pocitat jina 2
//    to uz mi prijde zbytecne slozity
// Qvadraticke prvky
// uzly co nejsou ve vertexech (jsou na hranach) se pocitaji jen kdyz jsou uvnitr patche, vede to na slozitejsi algoritmus, viz napr. oofem
void MEER :: MEER_SPR_nodal_recovery (void)
{
  long i, j, regid;
  NODE *node;
  
  GPA<const NODE> bp_recovered_nodes;
  
  // refined_node_values allocation
  if (valtype_1) {
    refined_node_values_1 = new Dvctr**[nREGIONs];
    for (i=0; i<nREGIONs; i++) {
      refined_node_values_1[i] = new Dvctr*[nNODEs];
      for(j=0; j<nNODEs; j++) {
    refined_node_values_1[i][j] = NULL;
      }
    }
  }
  if (valtype_2) {
    refined_node_values_2 = new Dvctr**[nREGIONs];
    for (i=0; i<nREGIONs; i++) {
      refined_node_values_2[i] = new Dvctr*[nNODEs];
      for(j=0; j<nNODEs; j++) {
    refined_node_values_2[i][j] = NULL;
      }
    }
  }
  
  Lvctr cavalues_1(nNODEs); // count added values == pocet pricteni hodnot do promenne refined_node_values_1; tou se pak refined_node_values_1 musi podelit
  Lvctr cavalues_2(nNODEs); // count added values == pocet pricteni hodnot do promenne refined_node_values_1; tou se pak refined_node_values_1 musi podelit
  
  for (regid=0; regid<this->nREGIONs; regid++) {
    // continue pokud v regionu neni zadny prvek
    if (REGIONs[regid].nvals <= 0)
      continue;
    
    // detekce basic patch a jeji applicability
    this->MEER_SPR_basic_patch_detection(regid);
    
    // loop over nodes with bp_appli[i] == 1 - basic patch recovery
    for (i=0; i<nNODEs; i++) {
      node = NODEs+i;
      
      this->MEER_SPR_bp_recovered_nodes_detection (node, bp_recovered_nodes);
      
      if (bp_recovered_nodes.give_size() == 0) continue;
      
      if (valtype_1) this->MEER_SPR_patch_recovered_nodes_compute(regid, bp_recovered_nodes, valtype_1, refined_node_values_1[regid], cavalues_1);
      if (valtype_2) this->MEER_SPR_patch_recovered_nodes_compute(regid, bp_recovered_nodes, valtype_2, refined_node_values_2[regid], cavalues_2);
    }
    
    // interpolovane hodnoty v uzlech zprumeruje
    for (i=0; i<nNODEs; i++) {
      if (NODEs[i].bp_appli) {
    if (valtype_1 && refined_node_values_1[regid][i] == NULL)  _errorr2("values at node %ld not approximated", i+1);
    if (valtype_2 && refined_node_values_2[regid][i] == NULL)  _errorr2("values at node %ld not approximated", i+1);
    
    // prumerovani
    if (valtype_1 && cavalues_1[i] > 1)  refined_node_values_1[regid][i]->dvdby(cavalues_1[i]);
    if (valtype_2 && cavalues_2[i] > 1)  refined_node_values_2[regid][i]->dvdby(cavalues_2[i]);
    
      }
    }
  }
}

void MEER :: MEER_SPR_compute_normals_at_elements (void)
{
  VectoR v1, v2;
  ELEMENT *elem;
  
  // normaly na elementech
  for (long i=0; i<nELEMs; i++) {
    elem = ELEMs+i;
    
    // continue if computed already
    if (elem->normal)  continue;
    
    // allocate
    elem->normal = new VectoR;
    
    if (elem->nnodes == 3) {
      v1.beP2P(&NODEs[elem->nodes[1]].coords, &NODEs[elem->nodes[0]].coords);  // vector from node[1] to node[0]
      v2.beP2P(&NODEs[elem->nodes[1]].coords, &NODEs[elem->nodes[2]].coords);  // vector from node[1] to node[2]
      
      elem->normal->beVectProduct(&v1, &v2);
    }
    else if (elem->nnodes == 4) {
      v1.beP2P(&NODEs[elem->nodes[2]].coords, &NODEs[elem->nodes[0]].coords);  // vector from node[2] to node[0]
      v2.beP2P(&NODEs[elem->nodes[1]].coords, &NODEs[elem->nodes[3]].coords);  // vector from node[1] to node[3]
      
      elem->normal->beVectProduct(&v1, &v2);
    }
    else _errorr("Unsupported number of nodes at element");
  }
}


// detekce basic patch a jeji applicability pro kazdy uzel
void MEER :: MEER_SPR_basic_patch_detection (long regid)
{
  long i, j, k, ncoeff;
  long nsp;
  long nsuperelems, superelems[64];
  //const FElement *felem;
  NODE *node;
  const ELEMENT *patchelem;
  VectoR *normal;
  
  int snode;
  Dvctr snodes;
  bool boundarytest;
  
  switch (SPR_boundary) {
  case MEER_SPRbpt_basic:  boundarytest = false;  break;
  case MEER_SPRbpt_nobnd:  boundarytest = true;   break;
  default: errol;
  }
  
  if (boundarytest)
    snodes.resize_ignore_vals(nNODEs);
  
  ncoeff = this->MEER_SPR_give_number_coefficients((MEER_SPRpatchType)REGIONs[regid].sprtype);
  
  //* compute normals at elements
  if (this->normal_at_nodes == false)
    this->MEER_SPR_compute_normals_at_elements();
  
  // smyce pres uzly a detekce patche
  for (i=0; i<nNODEs; i++) {
    node = NODEs+i;
    nsp = 0;                                              // number sampling/integration points
    this->give_superelems_to_node(i, nsuperelems, superelems);  // elements superior to node[i]
    
    // nulovani, alokace
    node->bp_appli = 0;                    // vynuluju
    node->bp_elems.realloc(nsuperelems);   // alokuju max velikost == pocet superelems, ale nektery vypadnou anzto nelezi v regionu
    node->bp_elems.resize_ignore_vals(0);  // resize na 0, z predchoziho regionu
    
    for (j=0; j<nsuperelems; j++) {
      
      if (regid == ELEMs[superelems[j]].regid) {
    node->bp_elems.add(superelems[j]);
    
    nsp += ELEMs[superelems[j]].nIPs1;
      }
    }
    
    if (nsp) {
      if (nsp >= ncoeff)   node->bp_appli = 1;
      else                 node->bp_appli = 2;
    }
    
    // detekce uzlu na hranici/boundary
    if (boundarytest && nsp) {
      if (node->bp_elems.give_size() == 1) {
    int nnodes = (ELEMs+node->bp_elems[0])->nnodes;
    if      (nnodes == 3) node->boundary = 3;
    else if (nnodes == 4) node->boundary = 2;
    else errol;
      }
      else {
    snodes.zero();
    
    // loop over bp_elems
    for (j=0; j<node->bp_elems.give_size(); j++) {
      patchelem = ELEMs+node->bp_elems[j];
      
      for (k=0; k<patchelem->nnodes; k++)
        snodes[patchelem->nodes[k]]++;
    }
    
    snodes[i] = 0;
    
    // loop over bp_elems
    for (j=0; j<node->bp_elems.give_size(); j++) {
      patchelem = ELEMs+node->bp_elems[j];
      
      snode = 0;
      for (k=0; k<patchelem->nnodes; k++)
        snode = snode + snodes[patchelem->nodes[k]];
      
      if (snode != patchelem->nnodes-1+2) break;
    }
    
    if (j != node->bp_elems.give_size())
      node->boundary = 1;
      }
    }
  }
  
  
  //
  for (i=0; i<nNODEs; i++) {
    node = NODEs+i;
    
    // pocitani normaly
    if (this->normal_at_nodes == false) {
      if (node->normal == NULL) {
    node->normal = new VectoR*[nREGIONs];
    memset(node->normal, 0, nREGIONs*sizeof(VectoR*));
      }
      
      if (node->bp_elems.give_size()) {
    normal = new VectoR;
    
    // loop over bp_elems
    for (j=0; j<node->bp_elems.give_size(); j++)
      normal->add(ELEMs[node->bp_elems[j]].normal);
    
    node->normal[regid] = normal;
      }
    }
    
    // pocitani rotation matrix
    if (node->normal[regid] != NULL) {
      double limit = 0.001;
      normal = node->normal[regid];
      VectoR v1, v2;
      bool normalize = false;
      
      double xa = normal->x < 0.0 ? -normal->x : normal->x;
      double ya = normal->y < 0.0 ? -normal->y : normal->y;
      double za = normal->z < 0.0 ? -normal->z : normal->z;
      
      // X
      if       (xa >= ya && xa >= za) {
    if (normal->x == 0.0) _errorr("collapsed normal");
    
    // normal == (1,0,0)
    if (ya/xa < limit && za/xa < limit) {
      node->rotmat(0,0) = 0; node->rotmat(0,1) = 1; node->rotmat(0,2) = 0;
      node->rotmat(1,0) = 0; node->rotmat(1,1) = 0; node->rotmat(1,2) = 1;
      node->rotmat(2,0) = 1; node->rotmat(2,1) = 0; node->rotmat(2,2) = 0;
    }
    else {
      v1.x = 0.0;
      v1.y =  normal->z;
      v1.z = -normal->y;
      normalize = true;
        }
      }
      // Y
      else if  (ya >= za && ya >= xa) {
    if (normal->y == 0.0) _errorr("collapsed normal");
    
    // normal == (0,1,0)
    if (xa/ya < limit && za/ya < limit) {
      node->rotmat(0,0) = 0; node->rotmat(0,1) = 0; node->rotmat(0,2) = 1;
      node->rotmat(1,0) = 1; node->rotmat(1,1) = 0; node->rotmat(1,2) = 0;
      node->rotmat(2,0) = 0; node->rotmat(2,1) = 1; node->rotmat(2,2) = 0;
    }
    else {
      v1.x = -normal->z;
      v1.y = 0.0;
      v1.z =  normal->x;
      normalize = true;
        }
      }
      // Z
      else {
    if (normal->z == 0.0) _errorr("collapsed normal");
    
    // normal == (0,0,1)
    if (xa/za < limit && ya/za < limit) {
      delete node->normal[regid];
      node->normal[regid] = NULL;
      
      // docasne, pak to smaz
      node->rotmat(0,0) = 1; node->rotmat(0,1) = 0; node->rotmat(0,2) = 0;
      node->rotmat(1,0) = 0; node->rotmat(1,1) = 1; node->rotmat(1,2) = 0;
      node->rotmat(2,0) = 0; node->rotmat(2,1) = 0; node->rotmat(2,2) = 1;
    }
    else {
      v1.x =  normal->y;
      v1.y = -normal->x;
      v1.z = 0.0;
      normalize = true;
        }
      }
      
      if (normalize) {
    v2.beVectProduct( normal, &v1 );
    
    v1.normalize();
    v2.normalize();
    normal->normalize();
    
    node->rotmat.copy_row(0,&v1);
    node->rotmat.copy_row(1,&v2);
    node->rotmat.copy_row(2,normal);
      }
    }
  } // loop over nodes
}


void MEER :: MEER_SPR_bp_recovered_nodes_detection (const NODE *node, GPA<const NODE> &bp_recovered_nodes)
{
  long i, j, elemid;
  bool supernode = false;
  bool cornertria = false;
  const NODE *patchnode;
  
  bp_recovered_nodes.resize_zero();
  
  // patch nejde vubec sestavit
  if (node->bp_appli != 1)  return;
  
  bool boundarytest;
  
  switch (SPR_boundary) {
  case MEER_SPRbpt_basic:  boundarytest = false;  break;
  case MEER_SPRbpt_nobnd:  boundarytest = true;   break;
  default: errol;
  }
  
  // uzel kolem ktereho je patch vytvorena se bude pocitat vzdy
  bp_recovered_nodes.add(node);
  
  // loop over bp_elems
  for (i=0; i<node->bp_elems.give_size(); i++) {
    elemid = node->bp_elems[i];
    
    // loop over element nodes
    for (j=0; j<ELEMs[elemid].nnodes; j++) {
      patchnode = NODEs+ELEMs[elemid].nodes[j];
      
      // bp_appli == 2  - add nodes
      if (patchnode->bp_appli == 2)
    bp_recovered_nodes.add_unique(patchnode);
      // bp_appli == 1
      else if (boundarytest) {
    // add node at boundary
    if (patchnode->boundary) {
      bp_recovered_nodes.add_unique(patchnode);
      if (patchnode->boundary == 3) cornertria = true;
    }
    // z patchnode se da spocitat tento node
    else {
      // pokud je tento node na boundary, tak se vyuzije moznost spocitat se z patchnode
      if (node->boundary)  supernode = true;
    }
      }
    }
  }
  
  if (boundarytest && node->boundary && supernode && !cornertria) {
    bp_recovered_nodes.resize_zero();
  }
}

// na patchi z elementu "bp_elems" se spocitaji koeficienty polynomicke approximace zvolene veliciny
void MEER :: MEER_SPR_patch_recovered_nodes_compute (long regid, GPA<const NODE> &bp_recovered_nodes, MEER_IPValues valtype, Dvctr **values, Lvctr &cavalues)
{
  long e, i, j, k, n, elemid;
  Dvctr *rhs, pol;
  Dmtrx eqs;
  const double *IP_values;
  
  long sprtype   = REGIONs[regid].sprtype;
  long nvalcomps = REGIONs[regid].nvals;
  long ncoeff    = this->MEER_SPR_give_number_coefficients(sprtype);
  
  //* main node
  const NODE *mnode = bp_recovered_nodes[0];
  
  //* rotation
  const Dmtrx *rotmat = &mnode->rotmat;
  const PoinT *coords;
  PoinT mnode_loc_coords, ip_loc_coords, node_loc_coords;
  bool mnormal = mnode->normal[regid] ? true : false;
  if (mnormal) {
    mnode_loc_coords.beRotatedPoint(rotmat, &mnode->coords);
  }
  
  
  //*  eqs[ncoeff,ncoeff] * coeffs = rhs   - hledame "nvalcomps" sad koeficientu
  
  eqs.resize_ignore_vals(ncoeff, ncoeff);
  eqs.zero();
  
  rhs = new Dvctr[nvalcomps];
  for (i=0; i<nvalcomps; i++) {
    rhs[i].resize_ignore_vals(ncoeff);
    rhs[i].zero();
  }
  
  // loop over patch element
  for (e=0; e<mnode->bp_elems.give_size(); e++) {
    elemid = mnode->bp_elems[e];
    
    for (i=0; i<ELEMs[elemid].nIPs1; i++) {
      IP_values = this->give_IPs1_values(elemid, i, valtype);
      
      if (mnormal) {
    ip_loc_coords.beRotatedPoint(rotmat, ELEMs[elemid].IPs1_coords+i);
    ip_loc_coords.sub(&mnode_loc_coords);
    coords = &ip_loc_coords;
      }
      else
    coords = ELEMs[elemid].IPs1_coords+i;
      
      this->MEER_SPR_give_polynom(sprtype, coords, pol);
      
      // compute ip contribution
      for (j=0; j<ncoeff; j++) {
    for (k=0; k<nvalcomps; k++)
      rhs[k][j] += pol[j] * IP_values[k];
    
    for (k=0; k<ncoeff; k++)
      eqs(j,k) += pol[j] * pol[k];
      }
    }
  }
  
  Dvctr *coeffs = new Dvctr[nvalcomps];
  
  if (! eqs.GaussSolve (coeffs, rhs, nvalcomps))
    errol;
  
  delete [] rhs;
  
  // spocitaji se hodnoty approximace v bp_recovered_nodes
  Dvctr vals(nvalcomps);
  
  const NODE *node;
  for (n=0; n<bp_recovered_nodes.give_size(); n++) {
    node = bp_recovered_nodes[n];
    
    if (mnormal) {
      if (n) {
    node_loc_coords.beRotatedPoint(rotmat, &node->coords);
    node_loc_coords.sub(&mnode_loc_coords);
      }
      coords = &node_loc_coords;
    }
    else
      coords = &node->coords;
    
    this->MEER_SPR_give_polynom(sprtype, coords, pol);
    
    vals.zero();
    
    for (i=0; i<nvalcomps; i++)
      for (j=0; j<ncoeff; j++)
    vals[i] += pol[j] * coeffs[i][j];
    
    
    // add vals to values
    if (values[node->id] == NULL) {
      values[node->id] = new Dvctr(nvalcomps);
      values[node->id]->zero();
      cavalues[node->id] = 0;
    }
    
    values[node->id]->add(&vals);
    cavalues[node->id]++;
  }
  
  delete [] coeffs;
}

int MEER :: MEER_SPR_give_number_coefficients (long sprtype)
{
  switch ((MEER_SPRpatchType)sprtype) {
  case MEER_SPRPT_2d:  return 3;
  default: _errorr("Unknown SPR patch type");
  }
  return -1;
}

void MEER :: MEER_SPR_give_polynom (long sprtype, const PoinT *coords, Dvctr &pol)
{
  switch ((MEER_SPRpatchType)sprtype) {
  case MEER_SPRPT_2d:  
    pol.resize_ignore_vals(3);
    pol[0] = 1.0;
    pol[1] = coords[0].x;
    pol[2] = coords[0].y;
    break;
  default: _errorr("Unknown SPR patch type");
  }
}


double MEER :: eval_interpol_fces_at_IPs2 (long elementID, int ipID, double *answer) const
{
  MEER_ElemDisplInterpol edi;
  double ncoords[3];
  double w;
  
  edi = this->give_element_EDI (elementID);
  
  w = this->give_IPs2_coords_natural (elementID, ipID, ncoords);
  
  switch (edi) {
  case MEER_EDI_2d_3n_linear:
    w = 2 * w;
    answer[0] = ncoords[0];
    answer[1] = ncoords[1];
    answer[2] = ncoords[2];
    break;
  case MEER_EDI_2d_4n_bilinear_1:
    w = 0.25 * w;
    answer[0] = 0.25 * (1 + ncoords[0]) * (1 + ncoords[1]);   
    answer[1] = 0.25 * (1 - ncoords[0]) * (1 + ncoords[1]);   
    answer[2] = 0.25 * (1 - ncoords[0]) * (1 - ncoords[1]);   
    answer[3] = 0.25 * (1 + ncoords[0]) * (1 - ncoords[1]);   
    break;
  case MEER_EDI_2d_4n_bilinear_2:
    w = 0.25 * w;
    answer[0] = 0.25 * (1 - ncoords[0]) * (1 - ncoords[1]);   
    answer[1] = 0.25 * (1 + ncoords[0]) * (1 - ncoords[1]);   
    answer[2] = 0.25 * (1 + ncoords[0]) * (1 + ncoords[1]);   
    answer[3] = 0.25 * (1 - ncoords[0]) * (1 + ncoords[1]);   
    break;
  default: _errorr ("Unknown edi");
  }
  
  return w;
}

double MEER :: give_element_characteristic_size (long elementID)
{
  MEER_ElemDisplInterpol edi = this->give_element_EDI (elementID);
  double             size = this->give_element_area (elementID);
  
  switch (edi) {
  case MEER_EDI_2d_3n_linear:      size = sqrt(2.0 * size);  break;
  case MEER_EDI_2d_4n_bilinear_1:
  case MEER_EDI_2d_4n_bilinear_2:  size = sqrt(      size);  break;
  default: _errorr ("Unknown edi");
  }
  
  return size;
}

} // namespace meerspace