MFEM v2.0
lininteg.cpp
Go to the documentation of this file.
00001 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
00002 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
00003 // reserved. See file COPYRIGHT for details.
00004 //
00005 // This file is part of the MFEM library. For more information and source code
00006 // availability see http://mfem.googlecode.com.
00007 //
00008 // MFEM is free software; you can redistribute it and/or modify it under the
00009 // terms of the GNU Lesser General Public License (as published by the Free
00010 // Software Foundation) version 2.1 dated February 1999.
00011 
00012 
00013 #include <math.h>
00014 #include "fem.hpp"
00015 
00016 void LinearFormIntegrator::AssembleRHSElementVect(
00017    const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
00018 {
00019    mfem_error("LinearFormIntegrator::AssembleRHSElementVect(...)");
00020 }
00021 
00022 
00023 void DomainLFIntegrator::AssembleRHSElementVect(const FiniteElement &el,
00024                                                 ElementTransformation &Tr,
00025                                                 Vector &elvect)
00026 {
00027    int dof = el.GetDof();
00028 
00029    shape.SetSize(dof);       // vector of size dof
00030    elvect.SetSize(dof);
00031    elvect = 0;
00032 
00033    const IntegrationRule *ir;
00034    if (IntRule)
00035    {
00036       ir = IntRule;
00037    }
00038    else
00039    {
00040       // ir = &IntRules.Get(el.GetGeomType(),
00041       //                    oa * el.GetOrder() + ob + Tr.OrderW());
00042       ir = &IntRules.Get(el.GetGeomType(), oa * el.GetOrder() + ob);
00043    }
00044 
00045    for (int i = 0; i < ir->GetNPoints(); i++)
00046    {
00047       const IntegrationPoint &ip = ir->IntPoint(i);
00048 
00049       Tr.SetIntPoint (&ip);
00050       double val = Tr.Weight() * Q.Eval(Tr, ip);
00051 
00052       el.CalcShape(ip, shape);
00053 
00054       add(elvect, ip.weight * val, shape, elvect);
00055    }
00056 }
00057 
00058 void BoundaryLFIntegrator::AssembleRHSElementVect(
00059    const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
00060 {
00061    int dof = el.GetDof();
00062 
00063    shape.SetSize(dof);        // vector of size dof
00064    elvect.SetSize(dof);
00065    elvect = 0.0;
00066 
00067    int intorder = oa * el.GetOrder() + ob;  // <----------
00068    const IntegrationRule &ir = IntRules.Get(el.GetGeomType(), intorder);
00069 
00070    for (int i = 0; i < ir.GetNPoints(); i++)
00071    {
00072       const IntegrationPoint &ip = ir.IntPoint(i);
00073 
00074       Tr.SetIntPoint (&ip);
00075       double val = Tr.Weight() * Q.Eval(Tr, ip);
00076 
00077       el.CalcShape(ip, shape);
00078 
00079       add(elvect, ip.weight * val, shape, elvect);
00080    }
00081 }
00082 
00083 void VectorDomainLFIntegrator::AssembleRHSElementVect(
00084    const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
00085 {
00086    int vdim = Q.GetVDim();
00087    int dof  = el.GetDof();
00088 
00089    double val,cf;
00090 
00091    shape.SetSize(dof);       // vector of size dof
00092 
00093    elvect.SetSize(dof * vdim);
00094    elvect = 0.0;
00095 
00096    int intorder = el.GetOrder() + 1;
00097    const IntegrationRule &ir = IntRules.Get(el.GetGeomType(), intorder);
00098 
00099    for (int i = 0; i < ir.GetNPoints(); i++)
00100    {
00101       const IntegrationPoint &ip = ir.IntPoint(i);
00102 
00103       Tr.SetIntPoint (&ip);
00104       val = Tr.Weight();
00105 
00106       el.CalcShape(ip, shape);
00107       Q.Eval (Qvec, Tr, ip);
00108 
00109       for (int k = 0; k < vdim; k++)
00110       {
00111          cf = val * Qvec(k);
00112 
00113          for (int s = 0; s < dof; s++)
00114             elvect(dof*k+s) += ip.weight * cf * shape(s);
00115       }
00116    }
00117 }
00118 
00119 void VectorBoundaryLFIntegrator::AssembleRHSElementVect(
00120    const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
00121 {
00122    int vdim = Q.GetVDim();
00123    int dof  = el.GetDof();
00124 
00125    shape.SetSize(dof);
00126    vec.SetSize(vdim);
00127 
00128    elvect.SetSize(dof * vdim);
00129    elvect = 0.0;
00130 
00131    int intorder = el.GetOrder() + 1;
00132    const IntegrationRule &ir = IntRules.Get(el.GetGeomType(), intorder);
00133 
00134    for (int i = 0; i < ir.GetNPoints(); i++)
00135    {
00136       const IntegrationPoint &ip = ir.IntPoint(i);
00137 
00138       Q.Eval(vec, Tr, ip);
00139       Tr.SetIntPoint (&ip);
00140       vec *= Tr.Weight() * ip.weight;
00141       el.CalcShape(ip, shape);
00142       for (int k = 0; k < vdim; k++)
00143          for (int s = 0; s < dof; s++)
00144             elvect(dof*k+s) += vec(k) * shape(s);
00145    }
00146 }
00147 
00148 void VectorFEDomainLFIntegrator::AssembleRHSElementVect(
00149    const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
00150 {
00151    int dof = el.GetDof();
00152    int dim = el.GetDim();
00153 
00154    vshape.SetSize(dof,dim);
00155    vec.SetSize(dim);
00156 
00157    elvect.SetSize(dof);
00158    elvect = 0.0;
00159 
00160    // int intorder = 2*el.GetOrder() - 1; // <-- ok for O(h^{k+1}) conv. in L2
00161    int intorder = 2*el.GetOrder();
00162    const IntegrationRule &ir = IntRules.Get(el.GetGeomType(), intorder);
00163 
00164    for (int i = 0; i < ir.GetNPoints(); i++)
00165    {
00166       const IntegrationPoint &ip = ir.IntPoint(i);
00167 
00168       Tr.SetIntPoint (&ip);
00169       el.CalcVShape(Tr, vshape);
00170 
00171       QF.Eval (vec, Tr, ip);
00172       vec *= ip.weight * Tr.Weight();
00173 
00174       vshape.AddMult (vec, elvect);
00175    }
00176 }
00177 
00178 
00179 void VectorBoundaryFluxLFIntegrator::AssembleRHSElementVect(
00180    const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
00181 {
00182    int dim = el.GetDim()+1;
00183    int dof = el.GetDof();
00184 
00185    shape.SetSize (dof);
00186    nor.SetSize (dim);
00187    elvect.SetSize (dim*dof);
00188 
00189    const IntegrationRule *ir;
00190 
00191    if (!IntRule)
00192       ir = &IntRules.Get(el.GetGeomType(), el.GetOrder() + 1);
00193    else
00194       ir = IntRule;
00195 
00196    elvect = 0.0;
00197    for (int i = 0; i < ir->GetNPoints(); i++)
00198    {
00199       const IntegrationPoint &ip = ir->IntPoint(i);
00200       Tr.SetIntPoint (&ip);
00201       const DenseMatrix &Jac = Tr.Jacobian();
00202       if (dim == 2)
00203       {
00204          nor(0) =  Jac (1,0);
00205          nor(1) = -Jac (0,0);
00206       }
00207       else if (dim == 3)
00208       {
00209          nor(0) = Jac (1,0) * Jac (2,1) - Jac (2,0) * Jac (1,1);
00210          nor(1) = Jac (2,0) * Jac (0,1) - Jac (0,0) * Jac (2,1);
00211          nor(2) = Jac (0,0) * Jac (1,1) - Jac (1,0) * Jac (0,1);
00212       }
00213       el.CalcShape (ip, shape);
00214       nor *= Sign * ip.weight * F -> Eval (Tr, ip);
00215       for (int j = 0; j < dof; j++)
00216          for (int k = 0; k < dim; k++)
00217             elvect(dof*k+j) += nor(k) * shape(j);
00218    }
00219 }
00220 
00221 
00222 void VectorFEBoundaryFluxLFIntegrator::AssembleRHSElementVect(
00223    const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
00224 {
00225    int dof = el.GetDof();
00226 
00227    shape.SetSize(dof);
00228    elvect.SetSize(dof);
00229    elvect = 0.0;
00230 
00231    int intorder = 2*el.GetOrder();  // <----------
00232    const IntegrationRule &ir = IntRules.Get(el.GetGeomType(), intorder);
00233 
00234    for (int i = 0; i < ir.GetNPoints(); i++)
00235    {
00236       const IntegrationPoint &ip = ir.IntPoint(i);
00237 
00238       Tr.SetIntPoint (&ip);
00239       double val = ip.weight*F.Eval(Tr, ip);
00240 
00241       el.CalcShape(ip, shape);
00242 
00243       add(elvect, val, shape, elvect);
00244    }
00245 }
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines