C++ Mesh::CalcMetrics方法代码示例

void Compute_dQdBeta_CTSE(Real* dQdB, SolutionSpace<Real>& space, Int beta)
{
  RCmplx perturb(0.0, 1.0e-11);

  Mesh<Real>* rm = space.m;

  SolutionSpace<RCmplx> cspace(space);
  Mesh<RCmplx>* cm = cspace.m;
  PObj<RCmplx>* cp = cspace.p;
  EqnSet<RCmplx>* ceqnset = cspace.eqnset;

  std::vector<SolutionSpaceBase<RCmplx>*> cSolSpaces;
  cSolSpaces.push_back(&cspace);

  SolutionOrdering<RCmplx> operations;
  std::string name = cspace.name;
  operations.Insert("Iterate " + name);
  operations.Finalize(cSolSpaces);

  Int cnnode = cm->GetNumNodes();

  Real* dxdb = new Real[cnnode*3];

  Get_dXdB(space.param->path+space.param->spacename, dxdb, rm, beta);

  //perturb complex part by dxdb*perturb
  for(Int i = 0; i < cnnode; i++){
    for(Int j = 0; j < 3; j++){
      cm->xyz[i*3 + j] += dxdb[i*3 + j]*perturb;
    }
  }
  //update xyz coords
  cp->UpdateXYZ(cm->xyz);
  //calculate metrics with new grid positions
  cm->CalcMetrics();
  
  //create a temporary operations vector with just a single space iterator on it
  //run solver
  Solve(cSolSpaces, operations);

  for(Int i = 0; i < cnnode; i++){
    for(Int j = 0; j < ceqnset->neqn; j++){
      dQdB[i*ceqnset->neqn + j] = 
	imag(cspace.q[i*(ceqnset->neqn+ceqnset->nauxvars) + j]) / imag(perturb);
    }
  }
  
  delete [] dxdb;

  return;
}

Real Compute_dObjdX_dXdBeta(SolutionSpace<Real>& space, Int beta)
{
  RCmplx perturb (0.0, 1.0e-11);

  Real dObjdBeta;
  RCmplx Obj;

  Mesh<Real>* m = space.m;
  EqnSet<Real>* eqnset = space.eqnset;
  Param<Real>* param = space.param;
  BoundaryConditions<Real>* bc = space.bc;

  SolutionSpace<RCmplx> cspace(space);
  Mesh<RCmplx>* cm = cspace.m;
  PObj<RCmplx>* cp = cspace.p;
  EqnSet<RCmplx>* ceqnset = cspace.eqnset;

  Int nnode = m->GetNumNodes();
  Int cnnode = cm->GetNumNodes();

  Real* dxdb = new Real[nnode*3];

  Get_dXdB(space.param->path+space.param->spacename, dxdb, m, beta);

  //perturb complex part by dxdb*perturb
  for(Int i = 0; i < cnnode; i++){
    for(Int j = 0; j < 3; j++){
      cm->xyz[i*3 + j] += dxdb[i*3 + j]*perturb;
    }
  }
  //update xyz coords
  cp->UpdateXYZ(cm->xyz);
  //calculate metrics with new grid positions
  cm->CalcMetrics();

  ComputeSurfaceAreas(&cspace, 0);
  Obj = Compute_Obj_Function(cspace);

  dObjdBeta = imag(Obj) / imag(perturb);
  
  delete [] dxdb;

  return dObjdBeta;
}

void Compute_dQdBeta_FD(Real* dQdB, SolutionSpace<Real>& space, Int beta)
{
  Real perturb = 1.0e-8;
  
  Int i ,j;
  EqnSet<Real>* eqnset = space.eqnset;
  Param<Real>* param = space.param;
  Mesh<Real>* m = space.m;
  PObj<Real>* p = space.p;

  Int nnode = m->GetNumNodes();

  std::vector<SolutionSpaceBase<Real>*> solSpaces;
  solSpaces.push_back(&space);

  //create a temporary operations vector with just a single space iterator on it
  SolutionOrdering<Real> operations;
  std::string name = space.name;
  operations.Insert("Iterate " + name);
  operations.Finalize(solSpaces);

  //compute dq/db with finite difference
  Real* dQdB_up = new Real[nnode*eqnset->neqn];
  Real* dQdB_down = new Real[nnode*eqnset->neqn];
  Real* Qorig = new Real[nnode*eqnset->neqn];

  for(i = 0; i < nnode; i++){
    for(j = 0; j < eqnset->neqn; j++){
      Qorig[i*eqnset->neqn + j] = space.q[i*(eqnset->neqn+eqnset->nauxvars) + j];
    }
  }

  Real* dxdb = new Real[nnode*3];

  Get_dXdB(space.param->path+space.param->spacename, dxdb, m, beta);

  //perturb by dxdb*perturb UP
  for(Int i = 0; i < nnode; i++){
    for(Int j = 0; j < 3; j++){
      m->xyz[i*3 + j] += dxdb[i*3 + j]*perturb;
    }
  }
  //update xyz coords
  p->UpdateXYZ(m->xyz);
  //recalculate metrics
  m->CalcMetrics();

  //perturb any parameters as required for sensitivities
  Perturb_Param(beta, *param, 1);
  //reset any interesting field which might depend on perturbations
  space.RefreshForParam();

  Solve(solSpaces, operations);
  for(i = 0; i < nnode; i++){
    for(j = 0; j < eqnset->neqn; j++){
      dQdB_up[i*eqnset->neqn + j] = space.q[i*(eqnset->neqn+eqnset->nauxvars) + j];
    }
  }

  //perturb by dxdb*perturb DOWN
  for(Int i = 0; i < nnode; i++){
    for(Int j = 0; j < 3; j++){
      m->xyz[i*3 + j] -= 2.0*dxdb[i*3 + j]*perturb;
    }
  }
  //update xyz coords
  p->UpdateXYZ(m->xyz);
  //recalculate metrics
  m->CalcMetrics();

  //perturb any parameters as required for sensitivities
  //back to original and then down
  Perturb_Param(beta, *param, -1);
  Perturb_Param(beta, *param, -1);
  //reset any interesting field which might depend on perturbations
  space.RefreshForParam();

  Solve(solSpaces, operations);
  for(i = 0; i < nnode; i++){
    for(j = 0; j < eqnset->neqn; j++){
      dQdB_down[i*eqnset->neqn + j] = space.q[i*(eqnset->neqn+eqnset->nauxvars) + j];
    }
  }

  //compute differences
  for(i = 0; i < nnode; i++){
    for(j = 0; j < eqnset->neqn; j++){
      //up
      dQdB_up[i*eqnset->neqn + j] -= Qorig[i*eqnset->neqn + j];
      dQdB_up[i*eqnset->neqn + j] /= perturb;
      //down
      dQdB_down[i*eqnset->neqn + j] -= Qorig[i*eqnset->neqn + j];
      dQdB_down[i*eqnset->neqn + j] /= -perturb;
    }
  }

  //average for central difference
  for(i = 0; i < nnode; i++){
    for(j = 0; j < eqnset->neqn; j++){
      dQdB[i*eqnset->neqn + j] = 
//.........这里部分代码省略.........

void Compute_dRdBeta_CTSE(Real* dRdB, SolutionSpace<Real>& space, Int beta)
{
  RCmplx perturb (0.0, 1.0e-11);

  Mesh<Real>* rm = space.m;

  Int nnode = rm->GetNumNodes();
  Int gnode = rm->GetNumParallelNodes();
  Int nbnode = rm->GetNumBoundaryNodes();

  std::ifstream fin;

  SolutionSpace<RCmplx> cspace(space);
  Mesh<RCmplx>* cm = cspace.m;
  PObj<RCmplx>* cp = cspace.p;
  EqnSet<RCmplx>* ceqnset = cspace.eqnset;
  Param<RCmplx>* param = cspace.param;
  
  Int cnnode = cm->GetNumNodes();
  Int cgnode = cm->GetNumParallelNodes();
  Int cnbnode = cm->GetNumBoundaryNodes();

  RCmplx* q = cspace.q;

  Int neqn = ceqnset->neqn;
  Int nvars = neqn + ceqnset->nauxvars;

  Real* dxdb = new Real[nnode*3];

  Get_dXdB(space.param->path + space.param->spacename, dxdb, rm, beta);

  //perturb complex part by dxdb*perturb
  for(Int i = 0; i < nnode; i++){
    for(Int j = 0; j < 3; j++){
      cm->xyz[i*3 + j] += dxdb[i*3 + j]*perturb;
    }
  }
  //update xyz coords
  cp->UpdateXYZ(cm->xyz);
  //calculate metrics with new grid positions
  cm->CalcMetrics();
  //compute gradient coefficients once
  ComputeNodeLSQCoefficients(&cspace);

  //perturb any parameters which we are getting sensitivities for
  Perturb_Param(beta, *cspace.param);
  //reset any interesting field which might depend on perturbations
  cspace.RefreshForParam();
  
  //the gaussian source is sometimes used to modify boundary velocities, etc.
  //pre-compute it for bc call
  if(cspace.param->gaussianSource){
    cspace.gaussian->ApplyToResidual();
  }

  //update boundary conditions
  cspace.p->UpdateGeneralVectors(q, nvars);
  UpdateBCs(&cspace);
  cspace.p->UpdateGeneralVectors(q, nvars);

  //now, compute gradients and limiters
  if(param->sorder > 1 || param->viscous){
    for(Int i = 0; i < (cnnode+cnbnode+cgnode); i++){
      ceqnset->NativeToExtrapolated(&q[i*nvars]);
    }
    cspace.grad->Compute();
    cspace.limiter->Compute(&cspace);
    for(Int i = 0; i < (cnnode+cnbnode+cgnode); i++){
      ceqnset->ExtrapolatedToNative(&q[i*nvars]);
    }
  }

  //compute spatial residual
  SpatialResidual(&cspace);
  ExtraSourceResidual(&cspace);
  for(Int i = 0; i < nnode; i++){
    for(Int j = 0; j < ceqnset->neqn; j++){
      dRdB[i*ceqnset->neqn + j] = imag(cspace.crs->b[i*ceqnset->neqn + j])/imag(perturb);
    }
  } 

  delete [] dxdb;
}

void ComputedQdBeta_HardsetBC(Int nnodes, Int* nodes, Real* dQdB, const SolutionSpace<Real>& space, Int beta)
{
  Real h = 1.0e-11;
  RCmplx I(0.0, 1.0);
  RCmplx perturb = h*I;

  Mesh<Real>* rm = space.m;

  SolutionSpace<RCmplx> cspace(space);
  Mesh<RCmplx>* cm = cspace.m;
  PObj<RCmplx>* cp = cspace.p;
  EqnSet<RCmplx>* ceqnset = cspace.eqnset;
  Param<RCmplx>* param = cspace.param;
  RCmplx* q = cspace.q;

  Int nnode = rm->GetNumNodes();

  Int neqn = ceqnset->neqn;
  Int nvars = neqn + ceqnset->nauxvars;

  Real* dxdb = new Real[nnode*3];

  Get_dXdB(space.param->path + space.param->spacename, dxdb, rm, beta);

  //perturb complex part by dxdb*perturb
  for(Int i = 0; i < nnode; i++){
    for(Int j = 0; j < 3; j++){
      cm->xyz[i*3 + j] += dxdb[i*3 + j]*perturb;
    }
  }
  //update xyz coords
  cp->UpdateXYZ(cm->xyz);
  //calculate metrics with new grid positions
  cm->CalcMetrics();
  //compute gradient coefficients once
  ComputeNodeLSQCoefficients(&cspace);

  //perturb any parameters which we are getting sensitivities for
  Perturb_Param(beta, *cspace.param);
  //reset any interesting field which might depend on perturbations
  cspace.RefreshForParam();
  
  //the gaussian source is sometimes used to modify boundary velocities, etc.
  //pre-compute it for bc call
  if(cspace.param->gaussianSource){
    cspace.gaussian->ApplyToResidual();
  }

  //update boundary conditions
  cspace.p->UpdateGeneralVectors(q, nvars);
  UpdateBCs(&cspace);
  cspace.p->UpdateGeneralVectors(q, nvars);

  for(Int i = 0; i < nnodes; i++){
    Int node = nodes[i];
    for(Int j = 0; j < neqn; j++){
      dQdB[i*neqn + j] = imag(q[node*nvars + j])/h;
      std::cout << "dqdb: " << i << ": -" << j << " " <<  dQdB[i*neqn + j] << std::endl;
    }
  }

}

Real Compute_dObjdBeta_CTSE(SolutionSpace<Real>& space, SolutionOrdering<Real>& operations, Int beta)
{
  Int i;  
  Real h = 1.0e-11;
  RCmplx I(0.0, 1.0);
  RCmplx perturb = h*I;

  Real dObjdBeta;
  RCmplx Obj;
  Mesh<Real>* rm = space.m;
  Param<Real>* param = space.param;
  BoundaryConditions<Real>* bc = space.bc;
  PObj<Real>* p = space.p;
  Int nnode = rm->GetNumNodes();

  SolutionSpace<RCmplx> cspace(space);
  Mesh<RCmplx>* cm = cspace.m;
  PObj<RCmplx>* cp = cspace.p;
  EqnSet<RCmplx>* ceqnset = cspace.eqnset;
  Int cnnode = cm->GetNumNodes();

  std::cout << "\n\nCOMPUTING dObj/dBeta_CTSE\n" << std::endl;

  Real* dxSurf = new Real[nnode*3];
  RCmplx* dxSurfC = new RCmplx[nnode*3];
  Compute_dXdB_Surface(space.param->path + space.param->spacename, rm, bc, dxSurf, beta);

  std::vector<SolutionSpaceBase<RCmplx>*> cSpaces;
  cSpaces.push_back(&cspace);

  for(i = 0; i < nnode*3; i++){
    //perturb complex part by displacement
    dxSurfC[i] = dxSurf[i]*perturb;
  }
  Int smoothingPasses = 1000;
  MoveMeshLinearElastic(cm, bc, dxSurfC, smoothingPasses);

  cm->CalcMetrics();

  //compute gradient coefficients once
  ComputeNodeLSQCoefficients(&cspace);

  SolutionOrdering<RCmplx> cOperations;
  cOperations = operations;
  cOperations.Finalize(cSpaces);

  //perturb parameters as required
  Perturb_Param(beta, *cspace.param);
  //reset any interesting field which might depend on perturbations
  cspace.RefreshForParam();

  //run solver
  Solve(cSpaces, cOperations);

  Obj = Compute_Obj_Function(cspace);
  dObjdBeta = imag(Obj) / h;
  
  std::cout << "\n\nComplex objective function value: " << Obj << std::endl;

  delete [] dxSurf;
  delete [] dxSurfC;

  return dObjdBeta;
}

void Compute_dRdBeta_FD(Real* dRdB, SolutionSpace<Real>& space, Int beta)
{
  Real perturb = 1.0e-8;

  EqnSet<Real>* eqnset = space.eqnset;
  Mesh<Real>* m = space.m;
  PObj<Real>* p = space.p;
  Param<Real>* param = space.param;
  Real* q = space.q;
  Int i;

  Int neqn = eqnset->neqn;
  Int nvars = neqn + eqnset->nauxvars;

  Int nnode = m->GetNumNodes();
  Int gnode = m->GetNumParallelNodes();
  Int nbnode = m->GetNumBoundaryNodes();

  Real* Qorig = new Real[nnode*nvars];
  Real* resOrig = new Real[nnode*neqn];
  Real* resUp = new Real[nnode*neqn];
  Real* resDown = new Real[nnode*neqn];

  //store original q
  for(i = 0; i < (nnode*nvars); i++){
    Qorig[i] = space.q[i];
  }

  UpdateBCs(&space);
  //compute spatial residual
  SpatialResidual(&space);
  ExtraSourceResidual(&space);
  for(Int i = 0; i < nnode; i++){
    for(Int j = 0; j < eqnset->neqn; j++){
      resOrig[i*eqnset->neqn + j] = space.crs->b[i*eqnset->neqn + j];
    }
  } 

  Real* dxdb = new Real[nnode*3];

  Get_dXdB(space.param->path+space.param->spacename, dxdb, m, beta);

  //perturb complex part by dxdb*perturb UP
  for(Int i = 0; i < nnode; i++){
    for(Int j = 0; j < 3; j++){
      m->xyz[i*3 + j] += dxdb[i*3 + j]*perturb;
    }
  }
  //update xyz coords
  p->UpdateXYZ(m->xyz);
  //calculate metrics with new grid positions
  m->CalcMetrics();
  //compute gradient coefficients once
  ComputeNodeLSQCoefficients(&space);

  //perturb parameters as required
  Perturb_Param(beta, *(space.param), 1);
  //reset any interesting field which might depend on perturbations
  space.RefreshForParam();

  //the gaussian source is sometimes used to modify boundary velocities, etc.
  //pre-compute it for bc call
  if(space.param->gaussianSource){
    space.gaussian->ApplyToResidual();
  }

  //update boundary conditions
  UpdateBCs(&space);
  //now, compute gradients and limiters
  if(param->sorder > 1 || param->viscous){
    for(Int i = 0; i < (nnode+nbnode+gnode); i++){
      eqnset->NativeToExtrapolated(&q[i*nvars]);
    }
    space.grad->Compute();
    space.limiter->Compute(&space);
    for(Int i = 0; i < (nnode+nbnode+gnode); i++){
      eqnset->ExtrapolatedToNative(&q[i*nvars]);
    }
  }

  //compute spatial residual
  SpatialResidual(&space);
  ExtraSourceResidual(&space);
  for(Int i = 0; i < nnode; i++){
    for(Int j = 0; j < eqnset->neqn; j++){
      resUp[i*eqnset->neqn + j] = (space.crs->b[i*eqnset->neqn + j]);
    }
  } 

  //perturb parameters as required, call twice (move back to original value then perturb down)
  Perturb_Param(beta, *(space.param), -1);
  Perturb_Param(beta, *(space.param), -1);
  //reset any interesting field which might depend on perturbations
  space.RefreshForParam();

  //the gaussian source is sometimes used to modify boundary velocities, etc.
  //pre-compute it for bc call
  if(space.param->gaussianSource){
    space.gaussian->ApplyToResidual();
  }
//.........这里部分代码省略.........