File: FieldFromAmplitudePhase.cpp

package info (click to toggle)
gmsh 4.8.4%2Bds2-3
  • links: PTS, VCS
  • area: main
  • in suites: bookworm
  • size: 87,812 kB
  • sloc: cpp: 378,014; ansic: 99,669; yacc: 7,216; python: 6,680; java: 3,486; lisp: 659; lex: 621; perl: 571; makefile: 470; sh: 440; xml: 415; javascript: 113; pascal: 35; modula3: 32
file content (148 lines) | stat: -rw-r--r-- 4,548 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
// Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/gmsh/gmsh/issues.

#include <stdlib.h>
#include "GmshGlobal.h"
#include "GmshConfig.h"
#include "GModel.h"
#include "MElement.h"
#include "OctreePost.h"
#include "FieldFromAmplitudePhase.h"

StringXNumber FieldFromAmplitudePhaseOptions_Number[] = {
  {GMSH_FULLRC, "Wavenumber", nullptr, 5.},
  {GMSH_FULLRC, "AmplitudeView", nullptr, 0.},
  {GMSH_FULLRC, "PhaseView", nullptr, 1.},
};

StringXString FieldFromAmplitudePhaseOptions_String[] = {
  {GMSH_FULLRC, "MeshFile", nullptr, "fine.msh"}};

extern "C" {
GMSH_Plugin *GMSH_RegisterFieldFromAmplitudePhasePlugin()
{
  return new GMSH_FieldFromAmplitudePhasePlugin();
}
}

std::string GMSH_FieldFromAmplitudePhasePlugin::getHelp() const
{
  return "Plugin(FieldFromAmplitudePhase) builds a complex field 'u' from "
         "amplitude 'a' (complex) and phase 'phi' given in two different "
         "'Views' "
         "u = a * exp(k*phi), with k the wavenumber. \n\n"
         "The result is to be interpolated in a sufficiently fine mesh: "
         "'MeshFile'. \n\n"
         "Plugin(FieldFromAmplitudePhase) generates one new view.";
}

int GMSH_FieldFromAmplitudePhasePlugin::getNbOptions() const
{
  return sizeof(FieldFromAmplitudePhaseOptions_Number) / sizeof(StringXNumber);
}

StringXNumber *GMSH_FieldFromAmplitudePhasePlugin::getOption(int iopt)
{
  return &FieldFromAmplitudePhaseOptions_Number[iopt];
}

int GMSH_FieldFromAmplitudePhasePlugin::getNbOptionsStr() const
{
  return sizeof(FieldFromAmplitudePhaseOptions_String) / sizeof(StringXString);
}

StringXString *GMSH_FieldFromAmplitudePhasePlugin::getOptionStr(int iopt)
{
  return &FieldFromAmplitudePhaseOptions_String[iopt];
}

PView *GMSH_FieldFromAmplitudePhasePlugin::execute(PView *v)
{
  double k = (double)FieldFromAmplitudePhaseOptions_Number[0].def;
  int aView = (int)FieldFromAmplitudePhaseOptions_Number[1].def;
  int phiView = (int)FieldFromAmplitudePhaseOptions_Number[2].def;
  std::string fileName = FieldFromAmplitudePhaseOptions_String[0].def;

  std::string name_model("");

  if(fileName == "") {
    Msg::Info("Could not find mesh file for interpolating U=A*exp(j*k*phi)."
              " Trying to use current model mesh, instead.");
    name_model = GModel::current()->getName();
    fileName = name_model + ".msh";
  }

  PView *va = getView(aView, v);
  if(!va) return v;
  PViewData *aData = va->getData();
  if(aData->getNumTimeSteps() != 2) {
    Msg::Error("Invalid number of time steps for AView, it must be complex!");
    return v;
  }

  PView *vphi = getView(phiView, v);
  if(!vphi) {
    Msg::Error("FieldFromAmplitudePhase plugin could not find PhiView %i",
               phiView);
    return v;
  }
  PViewData *phiData = vphi->getData();

  if(aData->hasMultipleMeshes() || phiData->hasMultipleMeshes()) {
    Msg::Error(
      "FieldFromAmplitudePhase plugin cannot be run on multi-mesh views");
    return v;
  }

  OctreePost *oA = nullptr, *oPhi = nullptr;
  oA = new OctreePost(va);
  oPhi = new OctreePost(vphi);

  GModel::current()->setVisibility(0);
  GModel *umodel = new GModel;
  umodel->readMSH(fileName);
  std::vector<GEntity *> _entities;
  umodel->getEntities(_entities);

  std::set<MVertex *> ve;
  std::map<int, std::vector<double> > dataR;
  std::map<int, std::vector<double> > dataI;

  for(std::size_t ent = 0; ent < _entities.size(); ent++)
    for(std::size_t ele = 0; ele < _entities[ent]->getNumMeshElements();
        ele++) {
      MElement *e = _entities[ent]->getMeshElement(ele);
      for(std::size_t nod = 0; nod < e->getNumVertices(); nod++)
        ve.insert(e->getVertex(nod));
    }

  for(auto it = ve.begin(); it != ve.end(); ++it) {
    double phi, ar, ai;
    std::vector<double> uR(1);
    std::vector<double> uI(1);
    oPhi->searchScalarWithTol((*it)->x(), (*it)->y(), (*it)->z(), &phi, 0);
    oA->searchScalarWithTol((*it)->x(), (*it)->y(), (*it)->z(), &ar, 0);
    oA->searchScalarWithTol((*it)->x(), (*it)->y(), (*it)->z(), &ai, 1);

    uR[0] = ar * cos(k * phi) - ai * sin(k * phi);
    uI[0] = ar * sin(k * phi) + ai * cos(k * phi);

    dataR[(*it)->getNum()] = uR;
    dataI[(*it)->getNum()] = uI;
  }

  delete oA;
  delete oPhi;

  PView *vu = new PView("FieldFromAPhi", "NodeData", umodel, dataR, 0.0, 1);
  vu->addStep(umodel, dataI, 1);

  if(name_model.empty())
    umodel->setName("fine");
  else
    umodel->setName(name_model);

  return vu;
}