Implementation of BoundaryIntegrator for ParMixedBilinearForm

[chikitkin]

Hi, everyone!

I want to implement second order absorbing boundary condition for electromagnetic modeling. In order to do that I have to define a ParMixedBilinearForm for the following expression (Jian-Ming Jin, The Finite Element Method in Electromagnetics, section 9.3 in the end):

$ \int_S \vec{V}_1 \cdot \nabla_t V_2 \, dS$ where $S$ - is a part of mesh boundary, $\vec{V}_1$ is a vector Nedelec function living on the whole mesh (H(curl)), $V_2$ is a scalar node based basis function living on the boundary S (e.g. H1(S)). Subscript $t$ denotes the tangential part (w.r.t. the surface) of the gradient.

I don't now how to do the following things:

1. How to define ParFiniteElementSpace for the $V_2$. It should be something like this auto fec = new H1_FECollection(pars.order, pmesh->Dimension()); auto fespace= new ParFiniteElementSpace(boundary_mesh, fec); but how to construct the boundary_mesh from the pmesh? It is the surface mesh on the part of the boundary of the space mesh pmesh.

2. Which integrator should I use for $ \int_S \vec{V}_1 \cdot \nabla_t V_2 \, dS$? Or, at least, which one should I choose as a starting point for implementation of my custom integrator? Will MFEM ''understand'' how to compute such a form for elements living on 3D mesh and surface mesh?

Thank you in advance for any advice!)

[mlstowell]

Hello, @chikitkin,

Performing operations on sub-meshes such as your boundary mesh is a topic that we have been discussing but we don't have support for this yet. However, I think we can find a solution in this case although it won't be optimal.

First, I would suggest defining the H1 field V_2 on the entire mesh. This should be straight-forward.

Second, use the MixedVectorGradientIntegrator to compute $ \int_S \vec{V}_1 \cdot \nabla_t V_2 , dS$. The trick is to add this as a boundary integrator like this:

ParMixedBilinearForm g(h1_fespace, nd_fespace);
g.AddBdrIntegrator(new MixedVectorGradientIntegrator, bdr_attr);

Where bdr_attr is an Array<int> object containing the boundary attributes of the surface S.

This will compute a rectangular bilinear form spanning the entire volume mesh but it will only be non-zero on the boundary S. Furthermore, it will automatically pick out the tangential portion of the gradient because the Nedelec basis functions on the boundary are designed so that the integrals of their normal components will be zero.

If you have any other questions please don't hesitate to ask. I'll be interested in learning how this turns out. Please consider sharing your test results if you can.

Best wishes, Mark

[chikitkin]

Thank you, @mlstowell ! This solves half of the problems)

I am now thinking about how to further construct a linear system. If I use the suggested bilinear form for entire mesh, then I have excess DoF for H1 - all except surface DoF. The simple solution with large overhead is to construct a block matrix

with a block equal to identity. Then all excess DoF in the solution would be zero. But may be these DoF can be eliminated from the unknown vector? May be I can define these excess DoF as true DoF for essential boundary and just set them to zero immediately? Could you please give the cue, what way is more suitable?

If I succeed I will indeed share my results for a test problem.

[mlstowell]

Hi, @chikitkin, I recently had to do something very similar. The approach I chose was to set the interior DoFs to zero using an essential boundary condition. The trick is to properly invert the list of true DoFs on the boundary to determine the true DoFs everywhere else. If you start with an Array<int> called bc_bdr_marker of length number-of-bdr-attributes containing zeros and ones indicating which attributes are on the boundary in question the following code should compute the true DoFs not on those boundaries.

      Array<int> bc_h1_tdof;
      H1FESpace_->GetEssentialTrueDofs(bc_bdr_marker, bc_h1_tdof);

      Array<int> bc_h1_tdof_marker;
      H1FESpace_->ListToMarker(bc_h1_tdof, H1FESpace_->GetTrueVSize(),
                               bc_h1_tdof_marker, 1);

      // Invert marker
      for (int i=0; i<bc_h1_tdof_marker.Size(); i++)
      {
         bc_h1_tdof_marker[i] = 1 - bc_h1_tdof_marker[i];
      }

      H1FESpace_->MarkerToList(bc_h1_tdof_marker, non_bc_h1_tdofs);

The array non_bc_h1_tdofs can then be passed to member functions such as ParBilinearForm::FormSystemMatrix and ParMixedBilinearForm::FormRectangularSystemMatrix to eliminate the appropriate rows and columns and, in the square case, place ones on the diagonal.

Best wishes, Mark

[chikitkin]

Thank you, @mlstowell ! I think that's exactly what I need, I will try this approach soon and report the results.

[chikitkin]

Hi, Mark,

I bumped into another intermediate problem. I do not understand how to correctly make EliminateEssentialBC in the case of a block matrix and nonzero Dirichlet conditions. Let's assume that I have the following block system

Here matrices $A, B, C, D$ come from assembling bilinear forms. For simplicity I consider the problem with the same spaces for u and v. Such problem for example results from rewriting complex problem via real bilinear forms. $u_2$, $v_2$ corresponds to DoF for essential B.C., so they are known. Let's consider the first block row of this system: We can move the known parts to the r.h.s.: and amend matrices $A$ and $B$ The identity block in $\tilde{A}$ will set b.c. automatically. The same can be done for $C$ and $D$:we set $C_22 = 0$, $D_22 = I$.

As far as I understand, function ParBilinearForm::EliminateEssentialBC do the following:

I also know that I can set zero block instead of identity passing dpolicy = DIAG_ZERO to this function. But if I apply EliminateEssentialBC twice for the same $b$ vector, it will rewrite $b_2$ with the b.c. values from the last call. Is there a some standard way to correct only a part of $b$ with EliminateEssentialBC or any similar function?

Sorry for large formulas, I help this comment will help someone else.

Aleksandr Chikitkin

[mlstowell]

Hi, Aleksandr,

You might be very interested to look at ParSesquilinearForm::FormLinearSystem where we struggle with this problem for a complex linear system defined as a 2x2 block system of real matrices as you mentioned.

There are two main tricks used to make this work in the complex case. The first is to call FormLinearSystem once for each block of the linear system but to use a temporary right hand side vector, equal to zero, in two of those calls to avoid overwriting the results of the previous calls. The second trick is to accept that your modified linear system will have not only A_22=I and D_22=I but also B_22=I and C_22=I. This works in the complex case because we end up with:

/I -I\
\I  I/

which is invertible. In general, if all of these were positive identities, it would fail. In the general case I think you would need to use dpolicy = DIAG_ZERO for the off-diagonal blocks. I have not worked through the details but I think the first trick would be sufficient in the general case. Specifically, use a series of temporary vectors, set equal to zero, for the right hand sides and combine the results after the calls to FormLinearSystem.

Best wishes, Mark

[chikitkin]

Thank you, Mark, it works!

[chikitkin]

Hi, Mark,

When I implemented construction of a block matrix similar to ParSesquilinearForm::FormLinearSystem, I found that this function differs from serial version SesquilinearForm::FormLinearSystem. Namely, in the serial version DiagonalPolicy::DIAG_ZERO is set for off-diag block:

blfi->SetDiagonalPolicy(mfem::Matrix::DiagonalPolicy::DIAG_ZERO);  
b_0 = 0.0;
blfi->FormLinearSystem(ess_tdof_list, x_i, b_0, A_i, X_0, B_0, false);
B_r -= B_0; 

while in the parallel version DIAG_ZERO is not set, and diagonal elements for true essential DoF are set to zero explicitly later:

b_0 = 0.0;
pblfi->FormLinearSystem(ess_tdof_list, x_i, b_0, A_i, X_0, B_0, false);
B_r -= B_0;
...
// Modify offdiagonal blocks (imaginary parts of the matrix) to conform
// with standard essential BC treatment
if (A_i.Type() == Operator::Hypre_ParCSR)
{
  HypreParMatrix * Ah;
  A_i.Get(Ah);
  hypre_ParCSRMatrix *Aih = *Ah;
  for (int k = 0; k < n; k++)
    {
      const int j = ess_tdof_list[k];
      Aih->diag->data[Aih->diag->i[j]] = 0.0;
    }
 }

I thought that the both ways do the same. But when I tried to use the approach from the serial version in my parallel code, I obtained in the end singular block matrix. And most likely this singularity is connected with the block of A_I corresponding to true ess-l DoF. Could you please explain, what is the difference between these two variants?

Aleksandr Chikitkin

[mlstowell]

Hello, Aleksandr,

If I recall correctly the serial and parallel implementations differ because, at the time, MFEM did not fully support the various DiagonalPolicy options in parallel. I also seem to recall that the zeroing out of the imaginary part in the parallel implementation may have been added later to improve solver performance.

It also seems to me that the serial and parallel implementations should be producing equivalent matrices. One of my colleagues also worked closely with this code and I'd like to ask for his advice. @psocratis, am I missing something? Are we correct in thinking the serial and parallel implementations should agree?

Best wishes, Mark

[psocratis]

It also seems to me that the serial and parallel implementations should be producing equivalent matrices. One of my colleagues also worked closely with this code and I'd like to ask for his advice. @psocratis, am I missing something? Are we correct in thinking the serial and parallel implementations should agree?

Hi @mlstowell , @chikitkin , indeed the serial and parallel implementations should agree (ex22 and 25 confirm this too). Mark, you are right, the difference between the two implementations was due to the fact that in the parallel case the diagonal policy was handled by Hypre and was explicitly set to one. The manual modification on the diagonal of the HypreParMatrix was done specifically for this purpose, i.e., for the serial and parallel codes to produce identical matrices. We also wanted to avoid extra non-zero entries in the off-diagonal blocks (for example in the case of a purely real-valued problem).

Best,

Socratis

[chikitkin]

Thank you, @psocratis, @mlstowell. I understand the difference now, but I still haven't managed to reproduce the solution obtained using ParSesquilinear form with my own implementation, where I explicitly construct a block matrix and solve the block system. I must have missed something, but I have no idea where is the error. I tried to make my code as close to the interior of ParSesquilinearForm::FormLinearSystem as possible. Apparently, the error in the construction of the block system, could you please look at it?

I checked, that sesquilinear form A = A_re + i A_im, where A_re and A_im are real forms from my code, and complex linear form b = 0 + i b_im, where b_im is real linear form. Here is the my code with block matrix (I eliminated non-important parts):

...
   // Read serial and parallel meshes.
   auto mesh = new Mesh(pars.fileNameMesh, 1, 0);
   //mesh->RemoveUnusedVertices();
   auto pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
   pmesh->ReorientTetMesh();

   // Define a parallel finite element space on the parallel mesh.

   auto fec = new ND_FECollection(pars.order, pmesh->Dimension());
   auto fespace = new ParFiniteElementSpace(pmesh, fec);
   // Determine essential boundaries 
   // (boundaries with Dirichlet b.c. (PEC - Perfect Electric Conductor)).
   Array<int> essentialBoundary = GetEssentialBoundary(pmesh, attributeList);
   if (myid == 0) {
      cout << "Essential boundary index: " << essentialBoundary << endl;
   }

   // Determine essential true dofs.
   Array<int> essentialTrueDofs;
   fespace->GetEssentialTrueDofs(essentialBoundary, essentialTrueDofs);

   ParLinearForm b_im(fespace);
   // ... Add integrators
   b_im.Assemble();

   // Construct form A_re
   ParBilinearForm* A_re = new ParBilinearForm(fespace);

   // int_V rot(\phi) rot(E)
   A_re->AddDomainIntegrator(new CurlCurlIntegrator(oneCoeffExt));
   // int_V -k_0^2 mu_r eps_r \phi E
   A_re->AddDomainIntegrator(
      new VectorFEMassIntegrator(coeffReal)
   );

   A_re->Assemble();

   // Construct form A_im
   ParBilinearForm* A_im = new ParBilinearForm(fespace);

   // int_V k_0^2 mu_r eps_r tand <\phi, E>
   A_im->AddDomainIntegrator(
      new VectorFEMassIntegrator(coeffImag)
   );
   // int_S k_0 <phi, E>
   A_im->AddBoundaryIntegrator(
         new VectorFEMassIntegrator(coeffBdrImag), 
         externalBoundary
   );

   // int_S (-1/2 k_0) <(rot E)_n, (rot phi)_n>
   A_im->AddBoundaryIntegrator(
      new CurlCurlIntegrator(beta),
      externalBoundary
   );
   A_im->Assemble();

   // Impose boundary conditions
   ParGridFunction * u_re = new ParGridFunction(fespace);
   ParGridFunction * u_im = new ParGridFunction(fespace);

   u_re->ProjectBdrCoefficientTangent(zeroVectCoeff, 
                                       essentialBoundary);
   u_im->ProjectBdrCoefficientTangent(zeroVectCoeff, 
                                       essentialBoundary);

   // Vector for solution
   Vector  X;
   X.UseDevice(true);
   X.SetSize(2 * tvsize);
   X = 0.0;
   //
   Vector X_r; 
   X_r.MakeRef(X, 0, tvsize);

   Vector X_i;
   X_i.MakeRef(X, tvsize, tvsize);

   // Allocate temprorary vector
   Vector b_0;
   b_0.UseDevice(true);
   b_0.SetSize(vsize);
   b_0 = 0.0;

   // Temprorary vector
   Vector X_0;

   // RHS vector for block system
   Vector B;
   B.UseDevice(true);
   B.SetSize(2 * tvsize);
   B = 0.0;

   Vector B_r; Vector B_i;
   B_r.MakeRef(B, 0, tvsize);
   B_i.MakeRef(B, tvsize, tvsize);

   // Temprorary vector;
   Vector B_0;
   // Check for memory leak;
   OperatorHandle A_11;
   OperatorHandle A_12; 
   OperatorHandle A_21; 
   OperatorHandle A_22; 

   HypreParMatrix* A_11_h = new HypreParMatrix();
   HypreParMatrix* A_12_h = new HypreParMatrix();
   HypreParMatrix* A_21_h = new HypreParMatrix();
   HypreParMatrix* A_22_h = new HypreParMatrix();

   // Compute A_11
   b_0 = 0.;
   // A_re->SetDiagonalPolicy(mfem::Matrix::DiagonalPolicy::DIAG_ONE);
   A_re->FormLinearSystem(essentialTrueDofs, *u_re, b_0, A_11, X_0, B_0);
   X_r = X_0;
   B_r = B_0;

   //Compute A_22
   // A_re->SetDiagonalPolicy(mfem::Matrix::DiagonalPolicy::DIAG_ONE);
   A_re->FormLinearSystem(essentialTrueDofs, *u_im, b_im, A_22, X_0, B_0);
   X_i = X_0;
   B_i = B_0;

   // Compute A_12: 
   b_0 = 0.;
   // A_im->SetDiagonalPolicy(mfem::Matrix::DiagonalPolicy::DIAG_ZERO);
   A_im->FormLinearSystem(essentialTrueDofs, *u_im, b_0, A_12, X_0, B_0, false);
   B_r -= B_0; 
   // A_12 *= -1.;

   // Compute A_21
   // A_im->SetDiagonalPolicy(mfem::Matrix::DiagonalPolicy::DIAG_ZERO);
   b_0 = 0.;
   A_im->FormLinearSystem(essentialTrueDofs, *u_re, b_0, A_21, X_0, B_0, false);
   B_i += B_0; 

   A_11.Get(A_11_h);
   A_21.Get(A_21_h);
   A_12.Get(A_12_h);
   *(A_12_h) *= -1.;
   A_22.Get(A_22_h);

   const int n = essentialTrueDofs.Size();
   // TODO: this part is different from SesquilinearForm::FormLinearSystem
   // I think, it is crucial only for CUDA
   for (int k = 0; k < essentialTrueDofs.Size(); k++){
      const int j = essentialTrueDofs[k];
      B_r[j] = X_r[j];
      B_i[j] = X_i[j];
   }
   // TODO: create A_ij as OperatorHandle and then cast them to HypreParMatrix
   hypre_ParCSRMatrix *A_12_ = *A_12_h;
   for (int k = 0; k < essentialTrueDofs.Size(); k++)
   {
      const int j = essentialTrueDofs[k];
      A_12_->diag->data[A_12_->diag->i[j]] = 0.0;
   }

   hypre_ParCSRMatrix *A_21_ = *A_21_h;
   for (int k = 0; k < essentialTrueDofs.Size(); k++)
   {
      const int j = essentialTrueDofs[k];
      A_21_->diag->data[A_21_->diag->i[j]] = 0.0;
   }

   // Synchronize memory

   B_r.SyncAliasMemory(B);
   B_i.SyncAliasMemory(B);

   X_i.SyncAliasMemory(X);
   X_r.SyncAliasMemory(X);

   Array2D<HypreParMatrix*> hBlocks;
   hBlocks.SetSize(2, 2);

   hBlocks(0, 0) = A_11_h;
   hBlocks(0, 1) = A_12_h;
   hBlocks(1, 0) = A_21_h;
   hBlocks(1, 1) = A_22_h;
   // Create BlockMatrix
   HypreParMatrix *H = HypreParMatrixFromBlocks(hBlocks);
   SuperLURowLocMatrix superlu_H(*H);
   SuperLUSolver solver(MPI_COMM_WORLD);
   solver.SetPrintStatistics(true);
   solver.SetOperator(superlu_H);   
   solver.Mult(B, X); 

   X_i.SyncAliasMemory(X);
   X_r.SyncAliasMemory(X);

   // Distribute solution
   u_re->Distribute(X_r);
   u_im->Distribute(X_i);

This code produces solution which differs from the correct solution computed using the following code with ParSesquilinear form:

   // 2.x Set up convention. 
   ComplexOperator::Convention conv =
      (pars.hermConv)
         ? ComplexOperator::HERMITIAN : ComplexOperator::BLOCK_SYMMETRIC;

   // 4. Read serial and parallel meshes.
   auto mesh = new Mesh(pars.fileNameMesh, 1, 0);
   //mesh->RemoveUnusedVertices();
   auto pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
   pmesh->ReorientTetMesh();

   auto fec = new ND_FECollection(pars.order, pmesh->Dimension());
   auto fespace = new ParFiniteElementSpace(pmesh, fec);

   // 7.a Determine essential boundaries.
   Array<int> essentialBoundary = GetEssentialBoundary(pmesh, attributeList);

   // 7.b Determine essential true dofs.
   Array<int> essentialTrueDofs;
   fespace->GetEssentialTrueDofs(essentialBoundary, essentialTrueDofs);

   ParComplexLinearForm b(fespace, conv);
   // Add integrators
   b.Assemble();

   // 9. Define solution vector u.
   ParComplexGridFunction u(fespace);
   u.ProjectBdrCoefficientTangent(zeroVectCoeff, zeroVectCoeff, 
                                       essentialBoundary);

   ParSesquilinearForm* a = new ParSesquilinearForm(fespace, conv);
   if (pars.pa) {
      a->SetAssemblyLevel(AssemblyLevel::PARTIAL);
   }

   a->AddDomainIntegrator(new CurlCurlIntegrator(oneCoeffExt), 
                          NULL);

   a->AddDomainIntegrator(
      new VectorFEMassIntegrator(coeffReal),
      new VectorFEMassIntegrator(coeffImag)
   );

   a->AddBoundaryIntegrator(
         NULL, 
         new VectorFEMassIntegrator(coeffBdrImag), 
         externalBoundary
   );  

   a->AddBoundaryIntegrator(
      NULL,
      new CurlCurlIntegrator(beta),
      externalBoundary
   );

   a->Assemble();

   OperatorPtr Ah;
   Vector B, X;
   a->FormLinearSystem(essentialTrueDofs, u, b, Ah, X, B);
   ComplexHypreParMatrix * AZ = Ah.As<ComplexHypreParMatrix>();
   HypreParMatrix * AC = AZ->GetSystemMatrix();
   SuperLURowLocMatrix SA(*AC);
   SuperLUSolver superlu(MPI_COMM_WORLD);
   superlu.SetPrintStatistics(true);
   superlu.SetOperator(SA);
   superlu.Mult(B, X);
   a->RecoverFEMSolution(X, b, u);

Thank you in advance!

[psocratis]

Hi @chikitkin, I believe I tracked down the issue. It's coming from the scaling *(A_12_h) *= -1.. Specifically the scaling propagates to A_21_h, probably because they are both associated with the same BilinearForm. You can either copy A_12_h to a new HypreParMatrix and then do the scaling (but this seems too much) or you can simply avoid doing the scaling at that stage and do it when you compute the monolithic matrix H by passing scaling coefficients. i.e.,

   Array2D<double> blocks_coeff;
   blocks_coeff.SetSize(2, 2);
   blocks_coeff(0, 0) = 1.0;
   blocks_coeff(0, 1) = -1.0;
   blocks_coeff(1, 0) = 1.0;
   blocks_coeff(1, 1) = 1.0;
   HypreParMatrix *H = HypreParMatrixFromBlocks(hBlocks,&blocks_coeff);

This seems to have fixed it for me.

I hope this helps. Please let us know if you still have the issue.

Best, Socratis

[chikitkin]

Thank you very much, @psocratis ! Now it works, I get the same solution.

[M8kmyday]

Back up at the top of this thread, there is the comment, "First, I would suggest defining the H1 field V_2 on the entire mesh. This should be straight-forward."

V_2 in this case is the component of the vector electric field normal to the surface. Assuming that the normal vector, nHat, is known, then V_2=nHat \cdot E. Since the goal is to build a BilinearForm, E is not yet known. How can V_2 be found from E for inclusion in the MixedVectorGradientIntegrator?

[M8kmyday]

Saying the same thing a different way similar to https://mfem.org/bilininteg/, a bilinear form such as this could work:

domain = ND, range = ND, coef = vector coefficient, operator = (\grad (vector_lambda \cdot vector_u), vector_v), continuous operator = \grad (vector_lambda \cdot vector_u), dimension = 2D, 3D

such that the operation on the 3D domain can be restricted to the boundary, where the ND elements automatically take the tangential component of the gradient, which is the required boundary condition.

For a little more background, this is the boundary condition (along with a vector mass term) that results from assuming a propagating 2D solution at the boundary. The goal is to use the boundary as a port representing a transmission line or waveguide feeding the 3D space. The general boundary condition is shown in https://github.com/mfem/mfem/issues/3610#issuecomment-1513587698, but the discussion there shows that the boundary condition is not supportable due to the derivative normal to the surface. By assuming a propagating solution with e^(-j*beta*z), the derivative is replaced by a simple multiplication and the resulting boundary condition includes the tangential gradient of the normal component of the electric field. A bilinear form such as that described above or an equivalent implementation would enable the complete implementation of the boundary condition under the given assumptions.

[mlstowell]

Hi, @M8kmyday,

I'm not sure I understand what you're proposing. If you hope to compute the tangential gradient of the normal component of a 3D Nedelec field using the bilinear form (\grad (vector_lambda \cdot vector_u), vector_v) integrated on the 2D boundary of a 3D mesh (with vector_u and vector_v in H(Curl) and vector_lambda equal to the surface normal), I'm afraid this is not going to work. The trouble is that a bilinear form computed on a surface only multiplies the DoFs contained within that surface. For a Nedelec field the surface DoFs correspond to basis functions which are mainly tangent to the surface. The normal component of a Nedelec field is primarily due to DoFs inside the domain contained in elements touching the surface but those DoFs are not on the surface.

Best wishes, Mark

[M8kmyday]

Hi Mark,

Thanks for the background. Rather than propose something, I'll just ask a question.

I am working with the electromagnetic weak integral form

As discussed in #3610, the surface integral cannot be set up in MFEM. The surface integral can be transformed by assuming a longitudinal dependence of e^(j*beta*z) to get

Here, z is being used, but in general the normal to the surface is the propagation direction. The tangential term can be added to a mixed bilinear form as

pmblfIm->AddBoundaryIntegrator(new MixedVectorMassIntegrator(*betaConst),*border_attributes);

I find this to be working well. I would like to include the gradient term, too, since Ez exists at the boundary, and it varies across the surface, so a gradient exists. Is there a way that I can include the gradient term?

[mlstowell]

Hi, @M8kmyday,

This cannot be done as a boundary integral for the reasons I mentioned earlier. The most accurate way I can think of would be to compute n.E everywhere in the domain as an H1 field and then take the gradient of that on the boundary. In place of n you could use any vector function which happens to equal the surface normal on the surface (as long as it doesn't deviate too much from that normal within the first layer of elements near the boundaries of interest). Use this vector function, n, in the MixedDotProductIntegrator within a MixedBilinearForm which maps ND to H1 fields. Then Solve an H1 mass matrix to compute the H1 field which most closely matches n.E. Finally you can use the MixedVectorGradientIntegrator as a boundary integrator to compute the surface integral of grad_t E_z dot W_m.

This would lead to a block linear system which would solve for E and n.E simultaneously. I have no idea how ill conditioned this might be.

Best wishes, Mark

[M8kmyday]

Thanks, Mark.

I think I see what you are saying and how to implement this. I have in general N ports with their own normals, ni. I will need to set up N calculations of ni.E over the volume of the domain to get the normal components needed for use in the boundary integrator. Instead of one volume calculation, I will have 1+N volume calculations in the block linear system which will greatly increase memory requirements and run time. Am I picturing this correctly?

If I am understanding this correctly, then I would be very hesitant to take this route.

[mlstowell]

Hi, @M8kmyday,

Yes, that's one way to do it. The more compact way is to make n a function of position. When it's near port i have it equal ni. It's value away from the ports is not important. The only values that matter are within the first layer of elements touching each port. Of course, attempting to support an arbitrary arrangement of ports of different sizes and orientations will make such a function difficult to implement but the benefit is that you only need to solve for n.E once.

Best wishes, Mark

[M8kmyday]

Just wanted to report in that I have Mark's suggestion up and running, and it seems to be doing the job. I'm not seeing any issues with ill conditioning so far with the small test jobs I am running for now. The needed custom vector coefficient for calculating n.E once turns out to be fairly simple. Thanks for the tips!