ベクトル（？）の要素の取り出し方

[Tohoko Tajima]

初めまして。最近OpenFOAMを勉強し始め、つまずいている部分があり質問させていただきます。

ほかの方が作ったライブラリ（https://github.com/vitst/dissolFoam）を修正して、自分のやりたいことを実現しようとしています。今躓いている部分は、以下のコードの修正です。おそらく私のOpenFOAMへの型の理解が問題で、ケースの問題設定等はあまり重要でないと思いますので、先にコードと質問を載せさせていただき、最後にやりたいこととコードの内容について書かせていただきます。

質問

• 以下のコードの tvalueという変数の Z成分だけに計算した値を入れるにはどうすれば良いのか教えていただきたいです。
• また、こういった型ごとの操作方法をみなさんはどうやって調べるのかを教えていただきたいです。

コード

このコードの tvalue  の1, 2つめの要素に値を入れずに０のままにしたく(太字下線部分)、いろいろ試したのですがエラーが出てしまいます。

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

template<class Patch>

template<class Type>

tmp<Field<Type> > CoupledPatchInterpolation<Patch>::faceToPointInterpolate

(

    const Field<Type>& ff

) const

{

    // Check size of the given field

    if (ff.size() != patch_.size())

    {

        FatalErrorIn

        (

            "tmp<Field<Type> > CoupledPatchInterpolation::"

            "faceToPointInterpolate(const Field<Type> ff)"

        )   << "given field does not correspond to patch. Patch size: "

            << patch_.size() << " field size: " << ff.size()

            << abort(FatalError);

    }

    tmp<Field<Type> > tresult

    (

        new Field<Type>

        (

//            patch_.nPoints(), pTraits<Type>::zero

            patch_.nPoints(), Zero

        )

    );

    Field<Type>& result = tresult.ref();

    List<Type> pointValue;

    pointValue.setSize( patch_.nPoints() );

   

    const labelListList& pointFaces = patch_.pointFaces();

    const scalarListList& weights = faceToPointWeights();

    forAll(pointFaces, pointi)

    {

        const labelList& curFaces = pointFaces[pointi];

        const scalarList& w = weights[pointi];

       

        Type tvalue = pTraits<Type>::zero;

        forAll(curFaces, facei)

        {

          tvalue += w[facei]*ff[curFaces[facei]]; // Z成分だけに対してこの計算をしたい

        }

        pointValue[pointi] = tvalue;

    }

   

    syncTools::syncPointList( mesh_, patch_.meshPoints(), pointValue, plusEqOp<Type>(), pTraits<Type>::zero);

   

    // normalization

    const scalarList& sumw = faceToPointSumWeights();

    forAll(pointFaces, pointi)

    {

      result[pointi] = pointValue[pointi] / sumw[pointi];

    }

   

    return tresult;

}

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

コード内の ff  や、 tvalue  をInfo<<でプリントしたところ、

tvalue = (0.002561356039 0.3089866244 -1.038888792)

ff[curFaces[facei]] = (-3.583025302e-06 7.391078324e-05 -0.0001789800755)

のようになります。 ff  の次元は[patchに含まれる点の数] x 3 (x, y, z成分),  tvalue  には３つの値 (x, y, z成分)が入っています。ここで、 tvalue  のx, y成分を０のままにするために、

forAll(curFaces, facei) // facei:0-3

{

tvalue[2] += w[facei]*ff[curFaces[facei]][2]; // w: scalar.  ff.size() = 19040

}

のような感じで、この計算をｚ成分だけにさせるようにしたいのですが、この書き方だと次のエラーが出ます（同じエラーがffのほうにも出ます）。

 (省略)..../lnInclude/CoupledPatchInterpolation.C:243:17: error: subscripted value is neither array nor pointer

 243 | tvalue[2] += w[facei]*ff[curFaces[facei]][2]; // w: scalar.  ff.size() = 19040

 | ~~~~~~^

[2]をつけて３番目の要素を取り出そうとするやり方は ff  と tvalue  の型（ const Field<Type>& ffとType tvalue = pTraits<Type>::zero; ）には使えないことはわかるのですが、代わりに何をすればよいのかわかりません。どうすればZ成分にだけこの計算をすることができるのか教えていただきたいです。

また、こういった型ごとの操作方法をみなさんはどうやって調べるのかを教えていただきたいです。

OpenFOAM-v1912、Ubuntu 22.04.1を使用しています。

最後に、質問内容的にあまり重要ではないかと思いますが、やりたいこととコードの内容について書かせていただきます。

やりたいこと

粗さを持つ２枚の板の隙間を流れる酸が、２枚の板を溶かす様子をシミュレートしたいです。そのために、ほかの方が作ったライブラリ（https://github.com/vitst/dissolFoam）の反応速度やメッシュ、メッシュの移動方法を修正しています。

コードの内容

上記のコードは、２枚の板の溶解体積を計算した後に、新しいメッシュの点の場所を決めています。このコードの前に、メッシュ状の各面の移動速度を決め、次に点の移動速度を決める際に点が囲まれた４つの面から移動速度をWeighted averageで決めています。

上記のコードを改変しようとしている理由は、このシミュレーションの結果をさらに別のシミュレーションに利用したく、その次の目的のためにメッシュの点がX, Y方向に移動されると困るからです。（この方法が数学的に、収束性的にやって良い事なのかは微妙ですがとりあえず試してみたいです。）

説明がとても長くなってしまいましたが、ご教授いただければ幸いです。

[ttsy shmz]

ttsyです。

テンプレートで書かれてるものをvector 決め打ちというのも不安ですが、以下の質問だけに答えます。

> tvalue  のx, y成分を０のままにするために

tvalue.x()=0;

tvalue.y()=0;

でxとy成分はゼロにできませんか？

iPadから送信

2023/02/17 9:04、Tohoko Tajima <msy...@gmail.com>のメール:



--
このメールは Google グループのグループ「OpenFOAM」に登録しているユーザーに送られています。
このグループから退会し、グループからのメールの配信を停止するには openfoam+u...@googlegroups.com にメールを送信してください。
このディスカッションをウェブ上で閲覧するには https://groups.google.com/d/msgid/openfoam/12fe829f-c916-4920-b14f-276f0157b10an%40googlegroups.com にアクセスしてください。

[Tohoko Tajima]

ttsyさん、

返信いただきありがとうございます！

tvalue.x(), tvalue.y()も試したのですが、エラーになってしまいました。他にも自分が思いつく限りの書き方は試したのですが、未だにうまくいっていません。

頂いた返信についてなのですが、「テンプレートで書かれてるもの」とはどういうことでしょうか？ケースディレクトリで使う .templateファイルのことは分かるのですが、このファイルはどの部分がテンプレートで書かれてるのでしょうか？

また、自分で試行錯誤する中で次の追加の質問があります。

1. 色々と試した際のエラーメッセージを見ると、tvalueやff[curFaces[facei]]（前のメールの太字下線部分の変数）がFoam::Vector<double>であると言われる場合と、const doubleであると言われる場合があり、結局これらの変数の型が何なのか分からず混乱しています。これらの変数の型を知る方法はエラーメッセージを見る以外にありますか？

　例えば、前のメールの太字下線部分の1行前でaという変数にff[curFaces[facei]]の値を入れようとすると、次のようなエラーがでます。

aをconst doubleにした場合：

　error: cannot convert ‘const Foam::Vector<double>’ to ‘const double’ in initialization //ff[curFaces[facei]]は‘const Foam::Vector<double>’と言っている？
  　247 |           const double a = ff[curFaces[facei]];
　　      |                        ^

aをVector<double>にした場合：

　error: conversion from ‘const double’ to non-scalar type ‘Foam::Vector<double>’ requested //ff[curFaces[facei]]は‘const double’と言っている？
　  247 |           Vector<double> a = ff[curFaces[facei]];

なぜこのように、場合によってff[curFaces[facei]]が違う型だと言われてしまうのでしょうか？

2. 上の質問に続き、tvalueやff[curFaces[facei]]がconst doubleであると言われる件についてなのですが、const doubleがベクトルのように3つの値を（(0.0 0.0 0.0)のような感じで）持つことができるのでしょうか？tvalueやff[curFaces[facei]]はInfo <<でプリントすると

　　tvalue = (0.002561356039 0.3089866244 -1.038888792)

　　ff[curFaces[facei]] = (-3.583025302e-06 7.391078324e-05 -0.0001789800755)

と3つの値が入っているのですが、自分でOpenFOAMの資料を調べたり適当な変数を作って試す限りでは、const doubleはscalarのように一つの値しか持っていないと思うので、どうしてtvalueやff[curFaces[facei]]がconst doubleであるとエラーメッセージで言われるのかがわかりません。

3. tvalueは "Type tvalue = pTraits<Type>::zero;"と宣言されているのですが、pTraitsというclassがどう使うものなのでしょうか？ネットで出てくる文章を読んでも使い道がよくわからなかったので質問させて頂きたいです。また、pTraitsにzeroではなく(0.0, 0.0, 1.0)のような感じで、各成分に違う値を入れることはできるのでしょうか？もしそれができれば、自分がやりたいことができるのではと思っています。

また長文かつ細かい質問になってしまいましたが、ご教授いただけると幸いです。

2023年2月17日金曜日 22:08:01 UTC-6 ttsy shmz:

[ttsy shmz]

ttsyです

わたしの言ったテンプレートは

です。

if(isA<vector>(tvalue)) tvalue.x()=tvalue.y()=scalar(0); で問題ないと思いますが、これでだめなんのですか？

iPhoneから送信

2023/02/21 3:36、Tohoko Tajima <msy...@gmail.com>のメール:



このディスカッションをウェブ上で閲覧するには https://groups.google.com/d/msgid/openfoam/7600132f-70bc-4853-95bd-3d89062edd87n%40googlegroups.com にアクセスしてください。

[Tohoko Tajima]

ttsyさん、

すぐに返信いただきありがとうございます。

また、テンプレートについて理解することができました。ありがとうございます！

if(isA<vector>(tvalue)) tvalue.x()=tvalue.y()=scalar(0); 

を試したのですが、以下のように、doubleは.x()と.y()を持っていないということでコンパイルエラーになってしまいます。

../../lnInclude/CoupledPatchInterpolation.C: In instantiation of ‘Foam::tmp<Foam::Field<Type2> > Foam::CoupledPatchInterpolation<Patch>::faceToPointInterpolate(const Foam::Field<Type2>&) const [with Type = double; Patch = Foam::PrimitivePatch<Foam::face, Foam::SubList, const Foam::Field<Foam::Vector<double> >&>]’:
dissolMotionPointPatchVectorField.C:2364:68:   required from here
../../lnInclude/CoupledPatchInterpolation.C:249:43: error: request for member ‘x’ in ‘tvalue’, which is of non-class type ‘double’
  249 |           if (isA<vector>(tvalue)) tvalue.x() = tvalue.y() = scalar(0);
          |                                              ~~~~~~~^
../../lnInclude/CoupledPatchInterpolation.C:249:56: error: request for member ‘y’ in ‘tvalue’, which is of non-class type ‘double’
  249 |           if (isA<vector>(tvalue)) tvalue.x() = tvalue.y() = scalar(0);
          |                                                                  ~~~~~~~^

しかし、プリントアウトすると、tvalueは(x, y, z)のように３つの値を持っています。

- なぜプリントアウトの見た目はベクトルで、エラーメッセージではdoubleと出てくるのでしょうか？

- この変数が持つメンバ関数を一覧表示する方法などがあれば教えていただきたいです。

- 自分で色々試していたところ、tvalueと全く同じように変数宣言したものは、以下のようにtvalueとcomptMultiplyで掛け算することができたので、これを(0.0, 0.0, 1.0)みたいにできれば自分のやりたい事は達成できると思っています。(コードの太字下線部以外は以前と同じです)

          Type tvalue - pTraits<Type>::zero;

          Type ttvalue - pTraits<Type>::one; //追加

          forAll(curFaces, facei)

          {

                    tvalue += w[facei] * ff[curFaces[facei]];

                    tvalue = cumptMultiply(tvalue, ttvalue); // 期待通りの計算結果になる

          }

ただ、pTraitsという型はzeroかoneしか設定できないのでしょうか？pTraitsというものが何なのか分からず(ドキュメントは読んだのですが、私の知識不足であまり理解できませんでした)、この方法もうまく行っていません。pTraitsに(0.0, 0.0, 1.0)のように値を設定できる方法があれば教えていただきたいです。

何度も同じ質問をしてしまって申し訳ないのですが、ご教授いただけると幸いです。

Tohoko Tajima

2023年2月20日月曜日 19:51:48 UTC-6 ttsy shmz:

[ttsy shmz]

ttsyです。

なるほど、と言うことは値が3個出力されてるけどvectorでは無いと言うことですね。

そうなるとdissolMotionPointPatchVectorField.Cを見ないと何とも言えないかな。

iPadから送信

2023/02/22 1:48、Tohoko Tajima <msy...@gmail.com>のメール:

このディスカッションをウェブ上で閲覧するには https://groups.google.com/d/msgid/openfoam/ddbea52d-de6b-47f5-88d6-c7556ebac0b2n%40googlegroups.com にアクセスしてください。

[ttsy shmz]

ttsyです

ユニバーサル参照を用いて

label n=0;

for(auto&& z : tvalue){

if(n==2) continue;

z=0;n++;

}

もダメですかね？

iPhoneから送信

2023/02/22 1:48、Tohoko Tajima <msy...@gmail.com>のメール:

このディスカッションをウェブ上で閲覧するには https://groups.google.com/d/msgid/openfoam/ddbea52d-de6b-47f5-88d6-c7556ebac0b2n%40googlegroups.com にアクセスしてください。

[Tohoko Tajima]

ttsyさん、

ユニバーサル参照を用いる方法も、tvalueがポインタなのでできなかったのですが、tvalueの値を変更する部分のみ、テンプレートで書かずに型によって別の関数を書くことで解決できました。

テンプレートのことや試すべきコードを教えていただかなければ、絶対自分では解決できなかったので本当に助かりました。ありがとうございます。

このファイルがテンプレートで書かれているので、コンパイル時に一旦doubleだと仮定してチェックされて(?)、それでこのコンパイルエラーが出続けていたんだと思います（2つめの返信で教えていただいたtvalueの型を用いたif文がどうして正しく動作してくれないのかはわかりませんが...）。

ffは、この関数に入る直前まではVector<double>で、.x(), .y()が使えたので、エラーメッセージでdoubleだと言われるのはおかしいと思い、tvalueを編集する部分だけテンプレートを使う方法に気づくことができました。

親身に教えてくださり本当にありがとうございました！

以下に一応変更部分を載せておきます。

template<class Patch>
template<class Type> // template program
tmp<Field<Type> > CoupledPatchInterpolation<Patch>::faceToPointInterpolate // used in dissolCase
(
    const Field<Type>& ff //cellMotionU for solubleWall

        Type fvalue = pTraits<Type>::zero;

        forAll(curFaces, facei) // facei:0-3
        {
          tvalue += w[facei] * ff[curFaces[facei]]; // w: scalar. ff: vector [3 x 1], ff.size() = 19040
        }

        AssignZeros(tvalue); //この関数を足した

        pointValue[pointi] = tvalue;
            Info << "tvalue in main = " << tvalue << nl;
    }

    // Info << pTraits<Type>::dim <<nl; // 3
    // Info << pTraits<Type>::typeName <<nl; // vector
    Info << "Using CoupledPatchInterpolation.C - Tohoko 3" << nl;    

    syncTools::syncPointList( mesh_, patch_.meshPoints(), pointValue, plusEqOp<Type>(), pTraits<Type>::zero);
   
    // normalization
    const scalarList& sumw = faceToPointSumWeights();
    forAll(pointFaces, pointi)
    {

      result[pointi] = pointValue[pointi] / sumw[pointi]; // pointValue is vector
      // result[pointi].x() = 0.0; // doesn't work
    }

    return tresult;
}

template<class Patch> //Typeがvectorだった場合
void CoupledPatchInterpolation<Patch>::AssignZeros
(
    vector& tvalue
) const
{
    tvalue.x() = 0.0;
    tvalue.y() = 0.0;
    Info << "tvalue in AssignZero = " << tvalue << nl;
}

template<class Patch>
void CoupledPatchInterpolation<Patch>:: AssignZeros  //Typeがdoubleだった場合
(
    double& tvalue
) const
{
    Info << "AssignZeros for scalar" << nl;
}

Tohoko Tajima

2023年2月21日火曜日 15:33:34 UTC-6 ttsy shmz: